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4 
From  a  cluster  expansion  of  the  pair  correlation  function  for  He   gas, 

the  term  linear  in  the  density  is  calculated,  and  the  third  virial  co- 
efficient obtained,  assuming  a  three-body  additive  potential  formed  by 
the  superposition  of  Lennard- Jones  pair  potentials. 

The  correlation  term  and  third  virial  coefficient  are  separated  into 
contributions  due  to  Boltzmann  statistics  and  exchange  effects.   The 
appropriate  conf igurational  density  matrix  elements  are  expressed  as 
functional  (Wiener)  integrals  which  are  approximated  by  multiple  Riemann 
integrals.   Numerical  results  are  obtained  by  Monte  Carlo  sampling  for  the 
temperatures  1  K  and  1.7  K. 

The  Boltzmann  pair  correlation  contribution  is  shown  to  agree  well 
with  a  so-called  "superposition  approximation"  based  on  two-particle  cal- 
culations at  both  temperatures.   While  contributions  due  to  statistics  to 
the  pair  distribution  function  are  small  at  1.7.  K,  they  are  seen  to  be 
significant  for  the  third  virial  coefficient  due  to  a  large  cancellation. 
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1.   INTRODUCTION 

4 
The  more  abundant  isotope  of  helium,  He  ,  has  long  been  the  focus  of 

considerable  experimental  and  theoretical  study.   Since  the  substance  is  of 

low  molecular  weight,  it  exhibits  many  properties  which  can  be  explained 

only  through  the  formalism  of  quantum  statistical  mechanics.   Basic  to  a 

4 
correct  understanding  of  the  properties  of  He   in  any  of  its  phases  is  a 

knowledge  of  the  intermolecular  interaction  potential.   Since  helium  mole- 
cules are  monatomic,  the  intermolecular  pair  potential  is  expected  to  be 
spherically  symmetric,  which  simplifies  greatly  calculations  in  which  it  is 
required.   While  attempts  at  accurately  calculating  the  pair  interaction 
potential  from  first  principles  over  a  complete  range  of  distances  have  thus 
far  met  with  limited  success,   studies  of  the  equilibrium  properties  of  the 
gas  have  confirmed  the  utility  of  certain  "empirical  potentials."   Compari- 
sons between  the  theoretical  and  experimental  temperature  dependence  of  the 
second  virial  coefficient  have  led  to  useful  parameter  values  of  such 
empirical  potentials  as  the  Lennard  Jones  (6-12)  [3]  and  the  modified 
Buckingham  (exp-6)  potentials  [4].   Further  information  about  the  pair 
interaction  as  well  as  the  three-body  interaction  potential  presumably 
could  be  learned  through  studies  of  the  third  virial  coefficient. 

While  a  considerable  number  of  calculations  of  the  second  virial  co- 
efficient exist  corresponding  to  a  variety  of  interaction  potentials,  cal- 
culations of  the  third  virial  coefficient,  which  in  general  involve  a 
quantum  mechanical  solution  to  the  three-body  problem,  have  been  limited. 
One  approach,  the  so-called  binary  collision  expansion  of  Lee  and  Yang  [5], 
has  been  used  by  Pais  and  Uhlenbeck  [6]  for  a  hard-sphere  potential  and 


See  Keller  [1]  for  a  survey  of  the  approaches  to  this  problem.   A  more 
exhaustive  account  is  contained  in  the  work  of  Hirshfelder  et  al.  [2], 


by  Larsen  [7]  for  a  square-well  potential  with  finite  repulsive  core.   The 
method  suffers  because  of  the  extreme  complexity  involved  in  evaluating 
terms  for  even  the  simplest  potentials,  and  the  necessity  of  calculating 
separately  the  two-  and  three-body  bound  state  energies  when  they  exist. 
In  addition,  the  approach  is  expected  to  be  applicable  only  for  low  tempera- 
tures.  Recently,  Larsen  and  Mascheroni  [8]  have  shown  how  the  third  virial 
coefficient  could  be  expressed  in  terms  of  three-body  phase  shifts  analogous 
to  the  usual  approach  to  the  second  virial  coefficient.   Currently,  this 
approach  is  limited  to  Boltzmann  statistics  only  and  does  not  take  into 
account  bound  states. 

Another  approach,  that  of  Jordan  and  Fosdick  [9],  has  been  to  formulate 
the  two-  and  three-particle  density  matrices  in  terms  of  path-integrals  and 
to  evaluate  the  integrals  by  Monte  Carlo  techniques.   This  approach  has  the 
advantage  of  having  no  inherent  limitations  as  to  the  nature  of  the  potential 
which  can  be  used  or  the  temperatures  to  be  considered.   Furthermore,  it  can 
yield  the  first  density-dependent  term  contributing  to  the  pair  correlation 
function  as  well,  a  three-particle  function  of  more  general  applicability 
than  the  third  virial  coefficient. 

Jordan's  calculations  were  based  on  the  assumption  of  Boltzmann  statis- 
tics and  an  additive  three-body  potential,  the  pair  interaction  being  the 
Lennard- Jones  potential.   With  decreasing  temperature,  computations  became 
increasingly  laborous,  and  the  sampling  error  tended  to  increase.   The 
lowest  temperature  achieved  was  5  K.   An  estimate  of  an  upper  bound  of  the 
exchange  terms  indicated  that  they  could  be  neglected  for  these  temperatures. 
The  upper  bound  estimate  was  analogous  to  one  previously  made  in  connection 
with  the  second  virial  coefficient  [10].   It  was  later  shown  by  Lieb  [11] 
that  a  suppression  of  effects  due  to  statistics  occurs  in  the  case  of  the 


second  virial  coefficient  which  is  much  stronger  than  was  indicated  by  the 
upper  bound.   The  suppression  of  statistics  was  shown  to  be  intimately 
related  to  the  strongly  repulsive  core  of  the  potential  which  limits  the 
possible  exchanges  to  those  which  do  not  occur  along  a  straight  line.   The 
exchange  contributions  to  the  third  virial  coefficient  are  likewise  expected 
to  be  suppressed  to  some  extent,  however,  no  direct  proof  of  this  has  been 
shown. 

The  present  work  is  an  extension  of  the  path  integral  approach  of 
Jordan  and  Fosdick  to  the  low-temperature  range  in  which  the  exchange 
effects  are  not  expected  to  be  negligible.   The  pair  distribution  function 
and  the  third  virial  coefficient  are  expressed  as  a  sum  of  a  term  corres- 
ponding to  Boltzmann  statistics  and  exchange  terms.   In  achieving  path  in- 
tegral calculations  for  the  low  temperatures  required,  an  attempt  is  made 
at  improving  the  approximation  scheme  and  the  sampling  procedure. 


2.   EQUATION  OF  STATE  FOR  HE   GAS 


2.1  Assumptions,  Units,  and  Conventions 

As  a  model  for  our  calculations  we  consider  a  gas  of  N  spinless  Bosons 

3 
contained  in  a  box  of  side  L  and  volume  L  .   It  is  assumed  that  L  is  large 

enough  so  that  the  size  and  shape  of  the  container  need  not  enter  into  the 

problem;  and  in  the  calculations  we  can  take  the  limit  L  —  «>,   N  -•  »  in  such 

N 
a  way  that  — r  -*  n.   In  order  that  the  above  assumptions  be  valid,  it  will 

L 

be  sufficient  that  an  actual  physical  gas  be  confined  to  a  container  the 

dimensions  of  which  are  much  greater  than  both  the  thermal  wavelength, 

,    (2TTQti2\   1/2     .  .         _  .   .       ,  TT  4 

A.  =     ( — K —  I    ,  and  the  range  of  the  interaction  between  any  two  He 

V  "He  / 

atoms.   We  also  assume  the  gas  is  in  a  state  of  thermal  equilibrium  and  as 

such  can  be  described  by  the  canonical  density  operator: 

-(% 

e 

where 

N 

""'i  i=lPi  +  V~N)  (2'l,1) 

N 
r  =  (r,,...,r  )  is  the  3N- dimensional  position  operator,  p.  is  the 

associated  3-dimensional  momentum  operator  for  particle  i  (=1,2,...,N)  and 


N 


VN(r  )  is  the  interaction  potential  for  the  N  particles.   Eventually 


we 


N 
will  assume  V„(r  )  to  be  in  the  form  of  a  superposition  of  two-body  Lennard- 

Jones  potentials, 


where 


V(r. .)  =  4a 
ij 


&) "  -  ft) 


(2.1.2) 


and 


r .  .  =   r  .  -  r  . 


-16  ° 

a  =  14.04  X  10    erg   a  =  2.56  A 


This  assumption  is  not  needed  until  Chapter  4.   The  values  for  a  and  a  were 
chosen  so  that  a  comparison  could  be  made  with  calculations  of  Jordan  [9]. 
They  were  obtained  by  deBoer  and  Michaels  [3],  and  provided  good  agreement" 
with  experiment  for  the  second  virial  coefficient  [12]. 

In  all  following  equations,  we  will  assume  units  such  that 

4 
ft  (Planck's  constant/2r0  =  irv,  (mass  of  He  molecule)  =  a    =  1.   In  these 

1     18  5 
units,  B  =  t-^  =  ~z —  where  k  is  the  Stephan-Boltzmann  constant  and  T  is  the 
kT    T 

o  40.97 

temperature  of  the  gas  in  K,  y  =   4a|3  =  — - —  and  the  thermal  wavelength, 

X  =  /2tt{3  ,  is  3.41//T.   We  will  generally  be  interested  in  the  conf  igurational 

matrix  elements  of  the  density  operator, 

<r'  |e     |r  > 
where 

If  )  S   Illf2' ' • • '5n> 

=  |l,2,...,N> 

is  the  eigenvector  of  the  position  operators  r.,    i  =  1,...,N.   The 
vectors  are  normalized  so  that 

,Ni  NN  _   3N,  ,N  N, 


/   ,  IN  •   INv        JIN.   ,  IN    IN. 

(r   |r  )  =  6   (r   -r  ) 


Calculations  by  Boyd,  Larsen,  and  Kilpatrick  [12]  actually  use 
a/k  =  10.22°K  in  accordance  with  [3].   The  value  adopted  here  and  by  Jordan 
yields  a/k  =  10.17°K,  the  discrepancy  apparently  being  due  to  improved 
measurements  of  k. 


The  superscript,  N,  will  generally  be  omitted  for  convenience  expect  in 
cases  where  it  is  necessary  for  clarity.   Occasionally,  the  superscript, 
op,  will  be  provided  to  distinguish  the  operator   from  its  eigenvalue. 

In  thermodynamic  equilibrium,  the  quantum  statistical  average  of  any 
observable,  8,  is  given  by 


-r-r     -1       "^ 

TeT  =  ZNL  tr(9e   1N)  (2.1.3) 


-pv . 


where  Z  =  tr(e    )  is  the  partition  function  of  the  system.   In  the  case 

of  identical  particles,  the  trace  must  be  taken  over  a  complete  set  of 

4 
properly  symmetrized  eigenstates.   Since  He  atoms  are  spin-zero  Bosons, 

the  states  must  be  simultaneous  eigenstates  of  the  permutation  operator, 

P,  with  eigenvalue  +1.  Accordingly  we  define  the  symmetrized  configurational 

eigenstates: 

|l,...,N>   h  -L.  Z   P|1,...,N> 

/NT  p 

where  these  kets  are  normalized  so  that  the  associated  wave  functions  are 
normalized  to  unity.  In  terms  of  symmetrized  configurational  states,  the 
partition  function  is 

ZN  =  nTJ  dr3N{l,...,N|e  ^|l,...,N>+      (2.1.4) 

where    the  N!    is   necessary   so   that  equivalent   states   related  by  a   permutation 
of  configurational  values  will  be   counted   only  once.      In  terms   of   the   un- 
symmetrized   kets   we   have 

ZN  "hS   dr3Nz    (1 N|Pe"PHN|l M> 


■  ^7  J  d3N  V£N>  •  <2-^> 


X3NN,  J  N 


whe  re 


W  (rN)  =  X3N  £  (1 n|Pe   ^|  1 N>      (2.1.6) 


N-         p 


is  known  as  the  Slater  sum.   The  diagonal  contribution  is  the  Boltzmann- 

_Q  U 

Slater   sum,    lOr   )    =   X     (l,...,N|e  |l,...,N>.      In   the   classical    limit 


N  B     N  -PVN(rN) 

Aim     W   (r   )    =      Aim     W~(r   )    =  e 

B  -  0  B  -»  0 


The  "generic"  conf igurational  distribution  function  for  h  particles, 
n,  (r  ),  is  defined  as  the  probability  density  of  the  occurrence  of  any  h 
particles  in  the  element  dr   about  the  point  r  =  (r. ,r« , . . . ,r,  ) .   Thus 


n.  (rh)  =   (Z         n   63(r°P  -  r,  ))  (2.1.7) 

h  ~        {L}  k=i     ~\    A 


where  X   represents  the  sum  over  all  permutations  of  h  particles  A  J„,.  .  ,  J, 
Uk]  l   2      h 

N! 
taken  out  of  the  total  N,  there  are  ...  ,  ^  ,  such  permutations,  and  each 

term  in  the  sum  is  equal.   Hence  according  to  equations  (2.1.3)  and  (2.1.5), 

equation  (2.1.7)  can  be  written, 

n    (rh)    = ^r^-   f    d3r  d3r  W    (rN)    .  (2.1.8) 

h~  (N-h)!X3NZNJ  h+1  NN~ 

In  the  case  h  =  2,  n,  (r  )  =  n„(r,,r0)  is  the  pair  distribution 

h  ~      2  ~1  ~2         r 

function.   For  h  =  1  we  have  the  particle  density,  n(r) .   In  the  limit  of 
large  volume,  V,  and  number  of  particles,  N,  it  is  well  known  that  for 
interactions  of  finite  range,  n(r)  —  n,  and  n9(r.,r9)  -•  n„  (  |  r_-r..  |  )  . 


2.2  The  Modified  Cluster  Expansion 

The  generic  distribution  function  n(r  )  can  be  expressed  as  an 
expansion  in  powers  of  the  fugacity,  z  =  ep  ,  where  u,  is  the  chemical 
potential.   The  so-called  "modified  cluster  expansion"  due  to  deBoer  [13], 

an  extension  of  the  cluster  expansion  by  Kahn  and  Uhlenbeck.  [14],  depends 

N 
on  a  transformation  of  the  Slater  sum,  W.,(r  ),  to  the  U  functions  defined 


by: 


W(rh)  =  U(rh) 

W(r\r  )  =  U(rh,r„)  +  U(rh)U(r  ) 

W(rh,rK,rx)  =  U(r\rK,rx)  +  U(rh,rR)U(rx)  +  U(rx,rK)U(rh) 

+  U(rh,rx)U(rK)  +  U(rh)U(rK)U(rx)       (2.2.1) 

where  we  have  dropped  the  subscripts  on  the  W's  for  convenience. 

Note  that  the  U-functions  are  defined  by  the  number  of  arguments  and 
not  the  dimensionality  as  with  W  functions.   Thus,  U(r  ,r   .)  and 
U(r, ,r„ , . . . ,r,  , 1 )  are  different  functions,  e.g.  for  h  =  2: 

U(r  ,r3)  =  W(r1,r2,r3)  -  W(r1,r2)W(r3> 

U(r1}r2,r3)  =  W^.r^^)  -  W(r1,r2)W(r3>  -  w(f2 ,53) W(5l) 
-  w(51'53)w(f2)  +  2W(Il)W(52)W(53)  ' 


There  are  thus  two  types  of  U-functions  implicitly  defined  above.   Those 
which  do  not  have  the  r  argument  are  the  functions  of  Kahn  and  Uhlenbeck. 
Those  which  do  include  the  r  argument  are  useful  in  obtaining  the  h- 
particle  generic  configurational  distribution  function,  n,  (r  ),  since  they 


become  small  whenever  the  particles  are  too  distant  from  one  another  or 
from  the  group  of  h  particles  at  location  r   to  interact.   This  property 
is  due  to  the  fact  that  the  W  function  tends  to  factor  into  products  of 
functions  involving  clusters  of  interacting  particles.   The  "interaction" 
distance  depends  upon  the  thermal-wavelength,  X,  as  well  as  the  range  of 
the  two-body  potential,  r_ ,  and  thus  increases  with  lower  temperatures. 
Thus  the  expansion  is  expected  to  converge  most  rapidly  for  sufficiently 
high  temperatures  and  for  short-ranged  potentials. 

The  actual  derivation  of  the  cluster  expansion  can  be  found  in 
de  Boer's  article.   The  result  is:' 

nK(rh)  =  \3hzh-1     I      ib<h)(rh)zX  (2.2.2) 

h  1=1   l        - 


where 

1 

x:x 


hT'{l  )  =  ,,.3i-3  I  U(I  •5w-i'"--»Iwi-i) 


x  d3rh+r--d3rh+x-i  (2-2-3) 


is  the  "modified  cluster  integral."   Equation  (2.2.2)  assumes  N  is  large; 
in  the  limit  as  V  becomes  large  so  that 

V  »  r   and  V  »  X3 
and 

N/V  -  n 

the  cluster  integrals  do  not  depend  on  the  size  or  shape  of  the  vessel 
containing  the  gas,  and  the  volume,  V,  may  be  regarded  as  infinite.   Under 


The  discrepancy  in  the  factor  of  X    is  due  to  a  difference  of  definition 
of  W(r^)  these  being  dimensionless . 
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such  conditions  the  functions  b*   (r)  are  independent  of  r  so  that 
b^Cr)  -  b^  and 

n,(r)  -\  2  Xb,zX  =  n  (2.2.4) 

1  ~    X3  X=l   X 


where  the  interchange  of  the  limit  and  the  summation  is  justified  as  long 
as  the  density  is  small  enough  that  the  particles  are  in  the  gas  phase 
[15]. 

In  the  case  of  h  =  2  we  have  the  pair  distribution  function  which 
depends  only  on  the  relative  distance  between  particles,  r19  =  | r_  -  r1 |  , 

00 


n  (r  r)-.!,  Z     ib[2)  (r.  _)Zi+1  .  (2.2.5) 

2  „1  ~2    x6   x=1   SL  12 


The  fugacity,  z,  can  be  eliminated  between  these  two  equations  yielding 

CD 


where 


n2(r12)  =7^  Vr12)(\3n)£ 


Yl  =  °'  Y2  =  bl2)(r12)'  Y3  =  2(b22)(^12>  "  2b2b[2)(r12)) 

Y4  =  3b^2)(r12)  -  12b2b^2)(r12)  +  (4b2  -  6b3)b<2)  (r^)   etc. 


Rewriting  the  above,  we  have 

n2(r12)  =  n2(n^0)(r12)  +  n^Cr^n  +  ...)  ,      (2.2.6) 


where 


and 


4° >(r   )  -  "2(5l,r2) 


11 

/  I  \  o 

n2   (r12^  =  J  d  r3^W3(5l'52'53')  "  w2(ri'r2)(W2(r2'r3) 

+  Vf3'fl}  "  l)^    •  (2-2.7) 

The  first  term  n~   (r)  involves  only  two-particle  interactions  whereas  the 
second  term  no   (O  involves  two-  and  three-particle  interactions. 

2.3   Expansion  of  the  Equation  of  State 

The  equation  of  state  of  a  gas  may  be  obtained  via  the  grand  canonical 
partition  function, 

00 

S  =  E  z  Z  (2.3.1) 

N=0 

The  pressure  of  the  gas  is  then  given  by 

PV  =  kT  log  E  .  (2.3.2) 

We  can  easily  expand  log  S  in  powers  of  the  activity,  z, 

00  CO  CO 

N  N  1  N        2 

log  S  =    log(l  +     I      Z   ziN)    =     S      z\   -  i   (  2      ziNZ    )Z 

N=l  N=l  N=l 

+  3    (  S     z   Z   )      -... 
N=l 


so 


00  i, 

PV   =   kT     Z      q,z 

i-1      l 


where 


A.  X 

q2   =   Z2   "  1  Zl   =  ~^6  I   ^l^^Vfl'^    '   Wl(Il)Wl(52^ 


=  ^el  d3rid3r2u2(5i'l2) 
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and   it   can  be   shown   that   in  general  [16] 

.3X 


1   r  a3JL  ,T  <  l\ 


Thus,  for  finite-ranged  potential, 

..   ^       h* 
*im  rr~  —  ~~ r 

v  -  «  v   xJ 


and 


CO  A 

kT  ^  ,   X 


Aim  P  =  ^  2  b.z   .  (2.3.3) 


As  with  the  pair  distribution  function  expansion,  we  eliminate  z  with  the 
expansion  for  n  (2.2.4)  and  obtain  a  power  series  in  the  molecular  density 
n.   In  order  to  obtain  the  familiar  virial  expansion,  we  express  P  as  a 
power  series  in  V   : 

|^-  =  A  +  B/V  +  C/V2  +  ...  (2.3.4) 

Substituting  from  equation  (2.2.4) 

CO  . 

V"1  =  ~     S   Xb.z 
NX3  1-1   X 

and  comparing  the  resulting  expansion  with  (2.2.10)  we  find 

A  =  1,  B  =  -  NX3b2,  C  =  N2X6(4b2  -  2b3)  . 

The  coefficients  A,  B,  C  are  the  first,  second,  and  third  virial  coefficients 
respectively.   Higher  order  virial  coefficients  may  be  obtained  in  the  same 
fashion.   It  is  to  be  noted  that  the  second  virial  coefficient  involves  no 
more  than  two-particle  interactions  and  the  third  virial  coefficient  in- 
volves no  more  than  three-particle  interactions.   We  may  express  them  in 
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terms  of  the  distribution  function  expansion  terms: 

2 
N  N       2 

B  =  ^  aQ   C  =  -  —  (a1-a0)  (2.3.5) 

whe  re 

00 

aQ  =  4tt  J  r2dr(l-n^0)(r))  (2.3.6) 

0 

CO 

aL  =  4tt  J  r2dr  n^1}(r)  .  (2.3.7) 

0 

Thus,  knowledge  of  the  first  two  terms  in  the  density  expansion  of  the  pair 
distribution  function  is  sufficient  to  determine  the  third  virial  co- 
efficient. 

2.4   Separation  of  Direct  and  Exchange  Terms 

We  now  use  (2.1.6)  to  express  n„   (r)  and  n„   (r)  in  terms  of 
diagonal  and  off-diagonal  configurational  density  matrix  elements.   For 
N  =  2  and  N  =  3  (2.1.6)  becomes 


W2(fl'f2)  =  X  L<12|e   Z|l2)  +  <2l|e    Z|l2>] 


=  n^0)(r12)        (2.4.1) 


9       ~^H3i  ~^H3 


W3(r15r2,r3)  =  \*[<123|e     | 123>  +  <213|e   J|l23> 

-pH  -PH 

+  <32l|e     | 123)  +  <132|e     | 123> 

-pH  -PH 

+  <23l|e     | 123>  +  <312|e     |l23>].    (2.4.2) 


The  off -diagonal  elements  represent  exchange  terms.   In  W_  there  is 
only  one  (pair)  exchange.   In  W  there  are  three  pair  exchange  and  two 
cyclic  exchanges.   We  can  represent  them  diagrammatically: 
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12  12  12  12  12 

(213|  <32l|  (132|  (23l|  <312| 


It  will  be   useful   to  arrange   the   integral,   n^      (r) ,   as    follows: 

n^Cr)   =  n<»<r)   +  n^fr)   +  n<l)<«)   +  n<l)(r)  (2.4.3) 


where 


Pl  -2 


n^lrjj)  =  J  d3r3[X9(123|e"PH3|l23>  -  X6<12|e"   2|l2> 

X  (Xb(23|e   J|23)  +  \b<3l|e   Z|3l)  -  1)}  (2.4.4) 

n^1}(r12)  =  2  J  d3r3X9<23l|  '3|l23>  (2.4.5) 

n^1}(r12)  =  J  d3r3{X9<213|e   3|l23)  -  \6<2l|e   2|l2) 

X  (Xb<23|e   Z|23)  +  \°<3l|e   Z|31>  -1)}  (2.4.6) 

np1)(r12)  =  2  I  d^3U9(l32|e   3|123)  -  \6<32|e   2|23) 

X  (Xb<12|e   2|12)  +  Xb(2l|e   Z|l2>)}  (2.4.7) 


and  where  the  factor  2  in  n^,   and  n!,   accounts  for  terms  in  which  r,«-»  r0 
which  are  equal.   The  terms  n   ,  xC      ,  and  np   (=  np   +  n   )  involve 
3-particle  direct,  cyclic  exchange,  and  pair  exchange  contributions  to  the 
pair  distribution  function,  respectively;  hence  the  subscripts.   Similarly, 
the  corresponding  contribution  to  the  second  and  third  virial  coefficients 


(2.3.6-7)  may  be  broken  up  into  direct  and  exchange  terms 
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D  ,   C  ,   P 

a,  =  a,  +  a,  +  a^ 


(2.4.8) 


where 


D   r  .3    (1),  v 

ai  =  fdrnc  (r) 


(2.4.9) 


3_  _(D 


ai  =  J  d  r  nc  (r) 


(2.4.10) 


J-rd3r(4l)(r)  +  4l)(r» 

J      Pl       P2 


(2.4.11) 
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3.   FUNCTIONAL  INTEGRALS  AND  THE  DENSITY  MATRIX 


3.1  The  Definition  of  the  Wiener  Integral 

The  subject  of  integration  over  function  space  with  respect  to  Wiener 
measure  has  been  treated  in  both  rigorous  and  non-rigorous  ways  by  many 
authors  [17-22]  and  can  be  approached  either  in  the  context  of  functional 
analysis  or  probability  theory.   The  present  approach  which  draws  from  a 
paper  by  Kac  [20]  represents  a  comprise  between  the  abstract,  but 
illuminating,  probabilistic  definition  and  a  simpler,  but  less  flexible, 
analytic  definition. 

The  Wiener  integral,  E[F[x(*)]}>  can  be  regarded  as  an  average  or 

expectation  of  the  functional,  F[x(*)]>  over  the  sample  space, 

Q  =  [continuous  functions  x(T):Te[0,t]  satisfying  x(0)  =  0],  with  the 

probability  measure  for  subsets  of  functions  being  assigned  in  the  following 

way.   Let  0  <  t,  <  t0  <  ...  t  <  t.   The  subsets  I  =  [x(-r):an  <  x(tJ  < 
12        n  u      1      1 

b,,...,a  <  x(t  )  <  b  }  H  Q  (referred  to  as  "measurable  rectangles"  or 
1      n      n     nJ  6 

"quasi-intervals")  are  assigned  the  probability 

b.         b 
1         n 

P(I)  =  J   dxx  ...  J   dxn  P(x1|0;t1)P(x2|x1;t2-t1)... 

a         a 

1         n 


P(x  x  ,  ;t  -t   , )  , 
n1  n-1  n  n-1   ' 


(3.1.1) 


where 


P(x,  x.  ;t,  -t  .)  = 
k1  j  k  y 


2tt(t,-t.) 


exp 


(vxi)  r 

2<vVj 


(3.1.2) 


By  simple  integration,  the  conditional  probabilities  satisfy 

CO 

P(xk,x.;VT.)  =  J  dx.P(xk|x.;VT.)P(x.|x.;T.-T.)     (3.1.3) 
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if  t.  <  t.  <  t.  ,  which  is  sufficient  to  insure  that  the  probability  assign- 
ment  on  I  is  consistent  with  the  assignment  for  any  equivalent  set,  e.g. 

I"  =  [x(T):a  <  x(t,)  <  b   -  cc  <  X(T«)  <*,...*     <  x(t  ) 
1—11  i         n      n 


<  b  }  0  Q 
nJ 


whe  re 


T{  +   TltT2...Tn 


By  a  theorem  due  to  Kolmogoroff,  there  exists  a  unique  probability  measure 
which  is  defined  over  a  Borel  field  of  subsets  of  Q,  and  which  agrees  with 
the  above  assignment  for  measurable  rectangles.   It  can  be  shown  that  the 
functions  determined  by  this  probability  measure  are  continuous  (i.e.,  dis- 
continuous functions  have  probability  zero).   These  functions  are  called 
Brownian  or  Wiener  paths.   The  integral  of  arbitrary  functionals  can  then 
be  defined  in  a  rigorous  way  [23].   For  our  purposes  it  is  not  necessary 
to  discuss  the  formal  definition  of  the  Wiener  integral.   We  will  however 
make  use  of  some  of  its  properties. 
Consider  the  functional, 

1      n 

F  [x(-)]   =  exp[-  -     E     V(x(t.))],      t.    =  -  t  (3.1.4) 

nLJ  rLn.riiJ  l        n 

i=0 

where  V(x)  >  0  is  a  continuous  function  of  x.   According  to  the  definition 
of  the  probabilities  (3.1.1),  the  average  or  Wiener  integral  of  this 
functional  is  given  by: 

E{Fn[x(.)]}  -  J  dxx  ...  J  dxnP(x1|0;i)  ... 


P(xJXn-l;n)exp[-  n  ;=  V(X1^  '     (3-1>5) 
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Now  consider  the  limit  as  n  -♦  ».   Since  V(x(t))  is  a  continuous  function 
of  t  (almost  everywhere)  it  follows  that 

n 


Xim  £  Z  V^c(t^  =  J  V(x(t)) 

n  -»  co   i=0         q 


dT   a.e. 


Further,  since  F  [x(#)]  is  bounded  and  measurable  it  follows  by  the  Lebesque 
dominated  convergence  theorem  [23]  that 

Xim     E{Fn[x(-)]}    -   E{    Xim     F^xC')]} 

n  — »  co  n  ~*  co 

or 

t  co  co  n-1 

E{exp[-  J   V(x(T))dr]]    =     Xim       J    ...  J   dx    . ..dx       II 

0  n  "*  "  -co     -oo  i=0 

n 

x  p<xi+ilVn)exp[-  n  E  v(xi>] 

i=0 


(3.1.6) 


where  x„  =  0. 


Of  particular  interest  in  quantum  statistical  mechanics  is  the 

conditional  Wiener  integral  which  will  be  defined  as  the  expectation  over 

the  subset  of  Wiener  paths  satisfying  x(t)  =  R.   It  can  be  shown  that  this 

conditional  expectation  can  be  expressed: 

t 
E{exp[-  J  V(x(T))dT]|x(t)  =  R]  = 


0 


n-1 


n  P(x.,1|x.;-) 
_•  n         i+l '  in 


03     »  i=     i+ii  i'n'         n 

Xim  J  ...  J  dx1...dxn_1   "      -    exp[-  "  L  V(x  )] 

n  -»  co_ro    _«  nl  >    '  1=0 

(3.1.7) 

where  x  =  R,  xn  =  0  and  the  factor  P(x  0;t)    is  necessary  so  that  the 
n       0  n1  J 

measure  is  normalized  to  unity.   It  can  be  shown  that  (3.1.7)  holds  for  all 
V(x)  bounded  from  below  and  possessing  only  a  finite  number  of 
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discontinuities.   The  appropriate  conditional  measure  for  measurable 
rectangles , 

I   =   [x(T):a1  <  xCt^   <  b1,...,an_1  <  xO^)   <  b^} 


0  <  t,...t      ,   <  t, 
1  n-1 


is   given  by 


n-1 


P(D    = 


Vi  .n0  p(-i+il-i^i> 

-[     **VJ        dVlX        P(x|0;t) 


(3.1.8) 


n-1 


where  x_  =  0  and  x  =  R.   The  continuous  paths  are  referred  to  as  conditional 
On  r 

Brownian  paths.   It  can  be  shown  that  they  satisfy  the  following  inter- 
polation formula  [24]: 

X(T)  •  «<TI)(T"-T?_+TX(T")(T-T-)  ^(T-T'y-T)^2   ^^ 


where  t"  <  t  <  t '  and  §  is  a  random  variable,  normally  distributed  with 
mean  0  and  variance  1. 

The  above  definitions  have  an  obvious  generalization  to  the  space  of 
3N-dimensional  functions,  x(t).  In  equation  (3.1.7)  we  simply  replace  x. 
by  x.  etc.,  and  define  the  conditional  probabilities  as 


P(x.  Ix.  ;t  -t  .)  =  oxt/o 

~kLj   k   J      [2TT(T  -T.)]3N/2 


exp 


<„v_V 

2(Tk-T.) 


2   1 


(3.1.10) 


2  . 


where  (x.  -x.)   is  the  square  of  the  distance  between  x,  and  x.  in  3N- 
dimensional  Euclidean  space. 


I'j 


3.2  Relation  Between  the  Density  Matrix  and  the  Wiener  Integral 

We  will  now  establish  the  important  connection  between  the  con- 
figurational  density  matrix  and  the  conditional  Wiener  integral.   The 
following  heuristic  argument  elucidates  the  relationship.   Consider  the 
Boltzmann  density  matrix  element: 

<r'|e^H|r> 
where 

H  ■  .2,  t.  +  V?  s  Ve>  +  V-r)  - 

1=1  1 


and  we  are  allowing  for  different  masses.   Now  we  can  write 

-fli      n   -fig 

e-PH  =  (e  n  }n  m     n   e  n  (32l) 

1=1 

since  H  commutes  with  itself.   Then  introducing  a  complete  set  of  eigen- 
states  of  r   ,  |x.),  between  each  operator,  we  obtain 

n-1  "££ 

<r'|e"pH|r>   =    T   dxf    ...    f   dx3N.    "a      <x.+Jen|x.)  (3.2.2) 

~    '  '-  «J         1  J        n-1    .    „     ~i+l'  '~i 

i=0 


where  we  have  defined  x_  =  r  and  x  =  r' .   We  now  consider  the  matrix 

~0       ~  ~n 

elements: 

-as 

(*i+ile  n  l5i>  =  ^i+J6^-  I  (tn  +  V%>  • 


Since  the  operators  T  and  VN  do  not  in  general  commute,  we  cannot  simply 
express  the  exponential  of  the  sum  as  the  product  of  exponentials: 
e      N  e  P    N#   However,  we  can  make  an  expansion  in  terms  of  the 
commutators,  assuming  |v(r)l  <  <=  and  (3  <  »: 
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e*p£-  !  <vvi =  e 


■p/n  TN      -p/n  VN 


X    (1  +  |(f)2[TN,VN]   +  0(^))  (3.2.3) 


where  0(— r)  involves  higher  order  commutators.   Thus, 
n 


<*i+l'e 


•p/n  H 


|e      |x.)  = 


■p/n  T. 


N, 


(^i+lle        l^i)e 


■P/n  V  (x  ) 

N  ~L  (1  +  0(-^))  . 
n 


(3.2.4) 


The  matrix  element  of  the  kinetic  energy,  T^,  is  easily  evaluated,  and  we 
finally  obtain: 


-PHi 


3N 


3N 


n-1  N 


1=0  k=l 


fa) 


3/2 


^  2 

exp["  W   (*i+l,k  "  *i,k}  3 


X  exp 


t-   I  V*^}  ti   +  0(^)]  • 


(3.2.5) 


Now,  with  the  transformation 


5i,k  -  Jl.k  ^  +  i  +  n  <Ik  "  *k>  k=1.---'N       <3-2-6> 


where  x.  =  (x.   ,x   ,...x.   ),  equation  (3.2.5)  then  becomes 

»**  JL        ^  i  j  1   •-  1  j  t.        *-  I  j  *  ' 
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n-1  3N/2  ? 

X     II     {(§-)  exp[-  f   (x       -x  )Z]} 

j=0        *"  ~J        ~J 

ft   n"1  i  1 

X  exp[-  *     2     VN(/R   y     +  r  +  ±<r '  -r))]  (l+OC^)) 

i=0  n 


(3.2.7) 


where 


y .   j 

ml       "n 

and  the  factor  (1  +  0(— r))  arises  from  the  product  of  n  factors  involving 

1  n 

(1  +  0(~)). 

n 
We  now  consider  the  limit  of  equation  (3.2.7)  as  n  -•  <=>.   The  error 

term  0(— )  is  assumed  to  approach  zero,  and  comparison  with  equation  (3.1.7) 

shows: 

/  it  -PH.  \    ,-3N  3?  r  3/2    r   "He  ,  ,    N2n. 
<r'|e  p  |r>  =  X    n  [m^      exp[-  ^  (&&   ]} 


k=l 


1 

X  E{exp[-ft  J  V  </0  ?(t)  +  r  +  r(r '  -r))]  |t](1)  =  0} 
0 

(3.2.8) 


where 


5(t)  i    —  T^t),...,—  T^Ct)  J  and  T1.(t)  (i  ^  1,...,N) 
ml  ~         "N  ~ 

are  three-dimensional  conditional  Brownian  paths.   Equation  (3.2.8)  was 
proved  by  Kac  [20,25].   The  only  requirement  on  V  (r)  is  that  it  be  bounded 
from  below  and  possess  only  a  finite  number  of  discontinuities. 
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ma 


A  useful  intuitive  interpretation  of  the  conf igurational  density 

trix,  (r'|e  | r) ,  is  provided  by  the  Wiener  integral  formulation  (3.2.8) 
We  consider  the  N  particles  as  undergoing  motion  from  conf igurational  points 
r  to  r'  as  in  Figure  1  (N  =  4) ,  the  trajectories  of  the  motion  being 


r 


j 


2' 


Figure  1.   Conditional  Brownian  Trajectories 
of  Four  Particles 


conditional  Brownian  paths  in  three  dimensions.   The  density  matrix 

r' 

i  -BH,  "PV? 

(r'je    r)  is  then  proportional  to  the  average  of  e   -   over  all  such 


paths  where  V~  is  the  time-averaged  interaction  potential  over  a  particular 

path.   Thus  (r'|e  ^  |  r)  <*   (e  ^  £)    where  the  brackets  on  the  right  hand 
denote  the  average  over  conditio n    Brownian  paths.   This  alternative 
notation  for  the  Wiener  integra'  .';:1  be  of  use  later.   It  is  to  be  noted 
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that  the  degree  of  randomness  in  the  path  is  controlled  by  the  thermal 
wavelength,  X.   Hence,  it  is  reasonable  to  attribute  the  randomness  of  the 
paths  to  the  Quantum  Mechanical  uncertainty  of  position  of  the  particles. 
In  the  classical  limit  as  X  ~*  0,  the  random  paths  become  straight  lines  and 
degenerate  into  points  for  r  =  r1,  and  because  of  the  "hard  core"  of  the 
interactions,  the  exchange  terms  go  to  zero.   In  the  direct  terms: 

<e-e^  >  -  e-ev(~r). 

3.3  Approximation  to  the  Wiener  Integral 
3.3.1  The  Mixed  Integration  Formula 

The  Wiener  integral  in  equation  (3.2.8)  cannot  in  general  be 

2 
evaluated  exactly  (there  are  some  exceptions,  e.g.  for  V(r)  =  r  ).   There- 
fore an  approximation  will  be  used.  We  begin  by  using  (3.2.8)  to  express 
equation  (3.2.2)  in  terms  of  Wiener  integrals: 

"X  2 

N  ' —   (r'-r  )  1 

X"3N     n     [m?/2   e   2P     ~k  ~k     }E{exp[-3    \   VG/&§(t)  +   r  +  T(r'-r))dT] 

k=l  o  ~  ~  ~     ~ 

X    |T)(1)   =  0} 

—  2 


_  r  ,  3N       ,  .  3n  n: :L  »  ,  («\y      . 
-  J  dxi   •••  J  dvi  nn  .n_t  [i^]     e*pl-\ 


2g/n 


i    , — 
X  E{exp[-  |  J  V(y|  |(t)  +  x.  +  T(xi+1-xi))dT]|T)(l)  =  0} 


where  xn  =  r,  and  x  =  r'  which  are  3N-dimensional  vectors  and  dx.   is 
~u   ~      ~n   ~  i 

3       3 
shorthand  notation  for  dx.  -...dx..  ...   We  make  the  transformation 

i,l     1,N 

x.  ,  -  x.  .   /E—  +  r,  +  -  (r,'  -  r,  ) 
~i,k   ~i,k  J  m,    ~k   n  ~k   ~k 
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and   finally  obtain, 

1 
E{exp[-0  J   V(/P|(t)   +  r  +  t  (r '  -r))dT]  |  T](  1)    -   0} 
0 

n-1 

n    p(x      |x  -±) 

OXT  TXT  '        1  ~1+1'~1      n 

-  r  dX3N      r  dx3N  — — , 

J    axl      '••   J    aVl  P(0   0;1) 


n-1  1  §(t) 

X     n      E{exp[-£JVN(/p(~-+zi  +  T(zi+1-z.)) 
i=0  o 


+  r  +  (^i)(r'-r))dT]|THl)  =  0}  (3.3.1) 


where 


§(t)  s  (t=-  T|  (T),...,~t-  \(t)) 
/ml  ~L        >n  ~N 

z   ■  (— —  x   ,(t),...,  — -  x  „(t))    i  =  0,1,..., n 


and 


xn  ,  -  x  .  =  0     k  =  1,2,. ..,N  . 
~0,k   ~n,k 

Equation  (3.3.1)  is  referred  to  as  the  "mixed  integration"  formula.   It  has 

the  appearance  of  a  truncation  of  (3.1.7)  except  that  the  function  of  the 

points  X-.X.....X   is  replaced  by  a  product  of  Wiener  integrals  of  functionals 

over  paths  between  the  points.   The  variables  z.  are  referred  to  as  "path 

break  points."   The  presence  of  —  multiplying  the  random  function  §(t) 

/n 

suggests  that  as  n  becomes  large,  the  statistical  variation  of  the  path 
(over  which  V  is  "time"  averaged)  becomes  small.   Similarly  as  (3  -♦  0,  this 
path  approaches  a  straight  line.   The  idea  of  expanding  the  conditional 
Wiener  integral  in  a  functional  Taylor's  series  about  /(3/n  T^t)  =  0  was 
suggested  by  Fosdick  [26],  and  a  modification  of  this  idea  was  employed  by 

*  r 

Jordan  in  his  doctoral  dissertation  |_  9]  . 


"Jordan  expanded  the  log  of  the  Wiener  integral. 
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3.3.2  The  "Simpson's  Rule"  Approximation 

Another  approach  of  Fosdick  and  Jordan  [27],  which  is  an  extension  of 
work  done  by  Cameron  [28]  to  the  case  of  conditional  Wiener  integrals,  is 
the  so-called  "Simpson's  rule"  approximation.   We  consider,  for  simplicity, 
the  one-dimensional  Wiener  integral,  E[f[ ( • )] | x(l)  =  0} .  A  function 

P(u,t)  is  chosen  so  that 

1 
E[P[x(-)]|x(l)  =  0]  =  \  I   du  P[p(u,-)]  (3.3.2) 

J  -1 

where  P„[x(0]  is  any  third  degree  polynomial  functional  of  the  form: 

1  1  1 

P3[x(-)]  =  hQ(T)  +  J  x(T)h1(T)dT  +  J  I    x(T1)x(T2)h2(T)dT1dT2 

0  0  0 

1  1  1 

+  J  I  I  x(Tl)x(T2)x(T3)h3(T)dTldT2dT3  '  (3-3'3) 
0  0  0 

It  can  easily  be  shown  that  equation  (3.3.2)  is  equivalent  to  the  conditions: 
(letting  r^  <   t2  <  t~) 

1 
\  J  du  p(u,T;L)  =  E[x(t1)|x(1)  =  0}  =  0  (3.3.4a) 

-1 
1 


|j  du  p(u,Tl)p(u,T2)  =  E[x(t1)x(t2)|x(1)  =  0] 

-1 


=  t1(1-t2)   ^  <  t2       (3.3.4b) 


1 

J 

-l 


2  J  du  p(u,T1)p(u,T2)p(u,T3) 

=  E[x(t1)x(t2)x(t3)|x(1)  =  0]  =  0         (3.3.4c) 


n 


where   the  moments,   E[   II      x(t,)|x(1)    =   0] ,    are   readily  calculated  by  direct 

i=l    L 

application  of  the  probability  measure  (3.1.8)  or  by  the  recursion  relation 
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(3.1.9).   It  is  not  difficult  to  verify  that  the  choice: 


P(u,t)  =  <   1-T   T  >  u  (3.3.5) 


-T          T    <    U 

u  >   0 

1-T       T    >    U 

-P(-u,t) 

u  <  0 

satisfies  conditions  (3.3.4).  Assuming  that  F[x(*)]  has  a  functional 
Taylor's  expansion  and  that  terms  of  order  higher  than  three  are  small,  we 
might  then  approximate  the  Wiener  integral  by  a  Rieman  integral: 

1 

E{F[x(-)]|x(l)  =  0}  S  \   J  du  F[p(u,-)]  (3.3.6) 

-1 

where  p(u,-r)  is  given  by  equation  (3.3.5).   This  method  of  choosing  an 
appropriate  family  of  parameterized  functions  so  that  the  functional 
integral  may  be  replaced  by  an  ordinary  integral  for  polynomial  functionals 
of  degree  less  than  four  is  analogous  to  the  Simpson's  rule  for  ordinary 
integrals  in  which  the  integral  can  be  replaced  by  a  finite  sum  for  poly- 
nomial functions  of  degrees  less  than  four.   Hence  the  approximation 
(3.3.6)  is  referred  to  as  the  "Simpson's  rule"  for  Wiener  integrals. 

The  "Simpson's  rule"  can  easily  be  extended  to  the  case  of  functionals 
in  more  than  one  variable.   We  simply  apply  the  following  rule:  replace  each 
independent  path  variable  T|.  (t)  of  the  functional  in  the  Wiener  integral  by 
the  function  p(u. ,t)  (3.3.5),  and  instead  of  taking  the  average  over  paths, 
take  the  average  with  respect  to  the  independent  random  variables  u. , 
i  =  l,2,...,n  where  each  u.  is  distributed  uniformly  over  [-1,1]. 

Konheim  and  Miranker  [29]  and  Fosdick  and  Jordan  [27]  showed  how  to 
generalize  the  approximation  for  Wiener  integrals  so  that  higher  order 
rules  could  be  obtained.   By  taking  linear  combinations  of  the  functions 
(3.3.5), 
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m 

(u.t)  = 
m 


(u,t)  =  S  c  p(u  ,t)  ,  (3.3.7) 

i=l 


and  replacing  9(u,t)  for  the  path  and  integrating  as  before  over  u.,...,u  , 
the  approximation  would  be  exact  for  polynomials  of  degree  2m+l,  providing 


2 
the  c.  are  the  roots  of  the  polynomial 


m    m-1  ,     m-1       ,  ,  lNm  ]     n  ,0  _  Qs 

z-z    +TTZ    -  ...  +  (-1)  — 7  =  0  .    (3.3.8) 

z .  m. 


It  can  be  shown  that  the  application  of  the  above  approximation  to 
the  functional, 


exp[-  |  J  V(x(t))<1t] 


yields  an  error  of  the  form  (8/n)   K  where  K  is  a  rather  complicated 

m        m 

integral  expression  involving  the  term 

1 

exp[-  |  J  V(s9m(u,T))dT] 
0 

where  s  is  integrated  over  the  range  [0,1].   In  the  case  of  the  Lennard- 

Y/(4n) 
Jones  potential,  the  exponential  term  is  bounded  by  e      ,  and,  because 

of  the  strongly  repulsive  part  of  the  potential,  is  probably  reasonably 

small  providing  9  (u,t)  is  real.   Unfortunately,  for  m  >  1  the  roots  of 

(3.3.8)  lead  to  complex  values  of  c.  so  that  the  K  can  be  quite  large, 

especially  for  calculations  in  which  the  repulsive  part  of  the  potential 

is  important.  Actual  path  integral  calculations  carried  out  by  this  author 

have  verified  this  problem.   For  this  reason,  the  present  work  does  not 

employ  these  higher  order  approximations. 

Applying  the  "Simpson's  rule"  to  each  of  the  Wiener  integrals  on  the 

right  side  of  equation  (3.3.1)  yields: 
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where 


and 


and 


E{exp[-  |  J   VN(/3    (i-  |(T)    +  fi   +  T(z.+1-z.)) 


/n 


+   r  +    (I±L-)(r'-r))dT]|Tl(l)    =   0] 
du3N  t 

2  0  /n 


+  t(z.+1-z.))   +  r  +   (I±i)(r'-r))dT]  +  3  (3.3.9) 


p'(u .  ,t)    =    ( p(u   .  ,t),.  .  .  ,—  p(u      ,t))      i   =   l,...,n, 

v/m1  /nijj 


£("ij,T)    =    (P(uiji'T)'P(uij2'T)'p(uij3'T^        J    =    1"--»N 


z.    =    ( x. .,..., x..T)        i   =   l,2,...,n-l 

-1         /mi  ~U  Vmj,  ~lN 


z      =    z      =   0 
~o       ~n 


u..,    x..  i   =   l,...,n;    1   =    1, . . . ,N  are   3 -tuples 


It   can  be   shown   that   the   error,   S,    is   0(— r).      There   are   n  such  integrals   in 

n 
the  mixed   integration   formula    (3.3.1);   hence    the    total  error   in  calculating 

,    ,N,    "0HN,    Nv    .         .         ,        nfl  v 
(r      |e  |r    )    is   of   order  C (— r) . 

n 
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4.   THE  COMPUTATIONAL  APPROACH 

4.1  Transformation  of  the  Matrix  Elements  to  the  Center  of  Mass  System 

The  density  matrix  elements  used  in  equations  (2.4.4-7)  can  be 
expressed  via  equation  (3.2.8)  in  the  form  of  Wiener  integrals  which  can  be 
approximated  by  the  mixed  integration  -  "Simpson's  rule"  combination.   The 
matrix  elements  are  of  dimensionality  6  and  9.   Since  we  are  assuming  the 
gas  is  contained  in  a  vessel  large  enough  so  that  its  size  and  shape  will 
not  enter  into  the  problem,  it  is  possible  to  transform  to  a  center  of  mass 
coordinate  system  so  that  we  can  reduce  the  dimensionality  of  the  density 
matrix  elements  to  3  and  6.   We  consider  now  3-particle  conf igurational 
density  matrix  elements.   Accordingly  we  make  the  transformation  to  the  R, 
r,  Q  coordinate  system  defined  by: 

5  =  I(£l+l2+53)     !r  =  £l  +  Zl   +  £3 

I   =  ll   -   5l        !r  =  ¥ll~l\> 

S  =  13   -   2(£l+f2)   !Q  =  3  £3  "  ?<Zl**Z>    '        (4'1-1) 


With  this  choice  of  position  and  momentum  coordinates,  the  corresponding 

operators  satisfy  the  usual  commutation  relations,  and  the  Jacobian, 

B(R  ,r  Q  ) 

r  ,     — ^ — ^ — r  ,  is  unity  (for  each  q  component).   Hence,  the  choice 

<nr]q>r2q'r3q; 

I R,  r,  Q)  =  1 1  2  3)  yields  properly  normalized  eigenkets  of  the  operators, 

N 
R   ,  r   ,  Q  P.   Assuming  V  (r  )  =  Z     V(r. .),  the  Hamiltonian  becomes 

i<j    1J 

P2    P2    P2 

H3  =  2m"  +  2m~  +  2m~  +  V^J   +   V^S  "  1   5>  +  V<S  +  2  ?      ^l'2> 
R     r     Q 


1    A  2 

lR   "»  mr  "  2  and  mQ  "  3 


whe  re  mR  =  3 ,  m  =  -r   and  r 
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We   next   consider   the   eigenvalues   of   the   permuted   kets,    P|l   2   3).      We 
obtain   them  by  operating  on  p| 1   2  3)   by  R      ,    r      ,    Q  P.      For  example    let 
P|l    2  3)   =    | 2   3    1);    then 

R°P|2   3   1>  -£   (r °P  +   r °P  +  r°P) | 2   3    1) 

■  5   (r2  +  r3  +  rx)|2   3    l)   =   r| 2   3    l) 

r°P|2  3    1)   -    (r°P   -    r°P)  1 2   3    1>   -    (^-r^  J2  3    l) 

=    (Q   -  \  r)|2  3    1) 
Q°P|2   3    1)   =[r°P   -\    (r°P+r°P)]|2   3    1) 

[rl   "   2    (r2  +  r3}^2  3    l) 
•  C-lj  -|r]|2  3   1)    . 


Hence , 


|  2   3   1>   =  c|R,    Q  -  \  r,    -  \  Q  -  |  r)  (4.1.3) 


where  c  is  determined  by  the  normalization  condition  and  is  equal  to  unity 
since  | R,  r ,  Q)  =  | 1  2  3) .   Similar  manipulations  yield  the  other  corres- 
ponding permuted  kets. 

|l  3  2>  «  |R,  Q  +  j  r,  -  \   Q  +  |  r)  (4.1.4) 

| 2  1  3>  =  |R,  -  r,  Q>  .  (4.1.5) 

Since  the  center  of  mass  coordinate,  R,  does  not  appear  in  the  Hamiltonian, 
the  Hamiltonian  may  be  broken  up  into  two  commuting  parts: 

H3  "  HR  +  "s"1 
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where 


and 


4 


2  2 

P  P 

h"1   =  ■—-  +  t5-  +  V(r)   +  V(Q   -  -^  r)   +  V(Q  +fr)    , 
3  2m  2m^  ~  I  ~  ~        ^  ~ 

r  Q 

Thus, 

(1   2   3|e        3|l   2   3)   =   <R,r,Q|e       *  e        r      |R,r,Q> 

rel 
=    <R|e        R|R>(r,Q|e        3      |r,Q> 

-pHrel 
=   33/2X"3(r,Q|e     '    3      |r,Q>    .  (4.1.6) 

Similarly, 

rel 
<2   1   3|e        3|l   2   3>   =   3J//\'J(-r,Q|e        J      |r,Q>  (4.1.7) 

-(3H3 


(2  3   l|e  |    12  3) 


-pHrel 
=   33/V3(Q  -  \  r,    -|q  -  |  r|e        3      |r,Q>  (4.1.8) 


(1   3   2|e        3|l   2  3) 

-6Hrel 
=  33/V3<Q  +  \   r,  -  \   Q  +  |  r|e  '  3   |r,Q>  .   (4.1.9) 

Similarly,  the  two-particle  density  matrix  elements  become: 

rel 

-j3H_  _  ,„      _  "BH,, 

<l'2'|e        l\\   2)   =  2J/V3<r'|e        Z      |r>  (4.1.10) 


where 


H"1  =   P2  +  V(r)    .  (4.1.11) 

2  r 
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4.2   Multiple-Integral  Approximation  of  the  Density  Matrix 

4.2.1   Three -Dimensional  Matrix  Elements 

We  now  combine  equations  (3.2.8),  (3.3.1),  and  (3.3.9)  to  obtain 
the  approximating  integrals  for  the  three-dimensional  density  matrix 
element  (4.1.10): 

-pHrel  ~(r'-r)  1 

23/2\-3<r'|e        2      |r>   =  e        4p  E{exp[-P  J   V(/2p   T1(t) 

0 

+  r  +  T(r,-r))dT]|T)(l)    =   0}  (4.2.1) 

-(r'-r)2  3n 


J   d^xn  J  ^Texp["   n  tZ=l 

X    f   V(/2£    (—  p(u     t)   +  x  +  T(x.-x.      )) 

0  yn 

+   r  +   (^i)(r'-r))dT]   +   O(^) 

n 


whe  re 


d^xn  =  [P(0|0;1)]"1     n     LP(xi+1|x.;  ^)]dx3.  .  .dx3^    ,        (4.2.2) 


and 


x     =   x     =   0 
~o       ~n 


P(U±,T)     S     (p(uil,T),p(ui2,T),p(ui3,T))      , 


and  we   have   expanded   the   product   of  n  3-dimensional   integrals   over   u  as 

one   3n-dimensional   integral   over   u.,u    ,...,u    ,    the   product   of  exponentials 

becoming  now   the   exponential   of  a   sum. 

du3n 
Since    '    du.        i    — - —  =    1,    the  multiple   integral    (4.2.1)    can  be   regarded 

as   an  average   of    the   exponential,    and   the   Monte   Carlo  method   can  be   employed 
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for  its  evaluation.   The  integral  over  t  can  be  evaluated  by  a  quadrature 
formula  or  by  an  exact  expression.   Since  p(u.,T)  is  a  discontinuous 
function  of  t,  the  point  of  discontinuity  for  the  q   component  being  u.  , 
it  will  first  be  necessary  to  break  up  the  integral  into  a  sum  of  integrals 
over  regions  in  which  p(u,t)  is  a  continuous  function  of  t.   To  do  this, 
we  consider  a  particular  value  of  u.  =  (u..,  u, „,  u.»).   Then,  defining 
the  time  partitions  t.  by:  Tn  =  0;  t.  =  |u..|,  j  =  1,...,3;  t,  =  1;  we  let 
k.  be  distinct  integers  e  {0,1,2,3,4}  which  satisfy:  k~  =  0,  k  =  4,  and 

0  S  Tk  <   Tk  <  Tk  <  \     <  Tk  H  l    ' 
0     1     2     3     4 

(We  need  not  consider  the  case  of  equalities  since  events  in  which  u's  are 

equal  occur  with  probability  zero.)   It  is  to  be  noted  that  t.  and  k.  both 

J      J 

depend  on  i;  we  will  not  explicitly  express  this  dependence,  however,  for 
notational  simplicity. 

The  integral  over  t  is  then  partitioned  in  the  following  manner: 

1  4   kJ+1    /55 

J  V(/2p  p(u,T)  +  ...)dT  =  2   J    V(  Jf-   p(u.,T)  +  ...)dT  . 

0  J=0  _ 

k. 


u.  . 

Let  S.  be  the  sign  function  of  u..,  i.e.  S.  =  i  1^>  ,  j  =  1,2,3.   Then, 
J  °  ij'       j    ju  I  »  J     »  « 

by  (3.3.5),  p(u   ,t)  =  p(S.t.,t)  =  S.p(t.,t).   Consider  now  the  interval 

(0,t   ].   Here  T  <  t.,  j  =  1,2,3;  so  from  (3.3.5),  p(t.,t)  =  -t;  hence 
Kl  J  J 

p(u. ,t)  =  (-S  t,-S  t,-S_t) .   Next  we  consider  the  interval,  (t,  ,t,  1: 
~  ~i         1    z    j  k  '  k„J 

we  have 

p(T   ,t)  =  1-t  p(u   ,t)  =  S   (1-t) 

Kl  kl       kl 

P(Tk  ,T)  =  -T  pd^  ,T)  =  -Sk  T 

P(T.   ,T)  =  -T  p(u   ,T)  =  -S   T 

R3  k3       k3 
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(i) 

U  1  l  V  ,       11      WC      UCLIUC      tllC       LU1CC      U  Lllli:  US  LUUd  i      vctiurs     (J 

j  =  0,1,2,3  by 

=  U'  Wj   'q  =  ^n~  xtL  V^'V5 

q  =  1,2,3       (4.2.3) 


and  so  on.   Finally,  if  we  define  the  three  dimensional  vectors  a 

~J 


°<T  =0'  toj^-v^^  V(q'k4);  j  =  lj2'3 


the    integral   over  t   in    (4.2.1)    takes    the    form: 

3  J'+1  t*\  /4X 

I        r  V(aU)    -   to'1'   +   ...)dT 

J-o  TJk>         -J 
J 


or 


1  3       kJ+1  r--^  (^ 

J   dT    ...    =     I       J  V(|uJlj;   +  U^\|)dT  (4.2.4a) 


0  J=0    T 


k 


where 


U,(ij)    =  a?i)   +  w.    ,    +  r  +  —   (r'-r) 
-1  ~j  ~i=l       ~         n       ~     ~ 

U^X)    =   -  a.(l)  +  A^.   +  -   (r'-r)  (4.2.4b; 

~z  ~j  ~i       n     ~     ~ 


and  U) 


.    =  /2g    x.  ,    AUJ.    -  U)     -   UU 
l  K    -l        ~i        ~i        ~i-l 


The    integral   over  t    is   now   in  a    form  where    it   can  be   evaluated  exactly 
(in   the   case   of   the    Lennard-Jones   potential)    or  by  means   of  a   quadrature 
method. 

4.2.2      Six-Dimensional  Matrix  Elements 

In   like  manner,   we   can  express   approximations    for   the   six-dimensional 

density  matrix  elements: 

,    -BH3el, 
(Q',r'|e  |Q,r>    . 


V, 


There  is  an  asymmetry  between  the  coordinates  Q  and  r  due  to  the  different 

2     1 
masses  -r  and  •=   associated  with  them;  accordingly,  the  associated  paths  will 

be  given  different  labels,  q(T)  and  T](t)  respectively,  and  the  corres- 
ponding "Simpson's  rule"  parameters  will  be  denoted  by  u.  and  v.  re- 
spectively.  It  can  be  seen  from  (3.2.8)  that  the  associated  Wiener  integral 
can  be  expressed  in  the  form: 

„3/2.6    r2  tt    ,   .2  ,  1  rr   .  ,   N2n 
3   \     exp[r-  — r  (Q'-Q)  +  "o  ~o  (r'-r)  ] 

*  x-    ~    ~  \ 

X  <Q',r' |e      |Q,  r)  (4.2.5) 

=  E{F("n(T),r,r,)F(Tl'(T),Q",Q,")F(Tl+(T),Q+,Q,+) 


where 


X  |T1(1)  =  q(l)  =  0} 


1 
F(|(t),x,x')  =  exp[-p  J  V(/2p  |(t)  +  x  +  T(x'-x))dT],   (4.2.6) 


+  +  + 

(|(t)  =  T)(t),T)'(t);   x  =  r,Q~;   x'  =  r',Q»") 

t      i 
T)  (t)  =  \   (q(T)  +  Tj(T»  ,  (4.2.7) 


t       1        t        1 

Q  -  Q  ±  2  f '    S   "  S  ±  2  5'  •  (4.2.8) 


It  will  be  noted  in  passing  that  the  Wiener  integral  associated  with  the 
3-dimensional  density  matrix  (4.2.1)  expressed  in  the  above  notation  is 

E{F(T](t),  r,  r')|Tl(l)  =  0}  . 


and 
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Returning  to  the  six- dimensional  case  and  applying  the  mixed  integration 
formula  (3.3.1)  and  the  "Simpson's  rule,"  we  obtain 

E{F(?,(T),r,r,)Fr+(T),Q+,Q,  +  )F(Tl"(T),Q",QO|T1(l) 

=  q(i>  =  0}  =  r  du.    r  d^    r  ^  r  ^1 

~       J   J   ^xn  J   ^yn  J   3n  J   3n 

R   n   - 
X  exp[-  *  E   (V(p(u.,T),x.,r,r') 

i=l 
+  V(p+(u.,y.,T),x^,Q+,Q1+) 

+  ^(£~(~i'~i'T)'*i'2~'3'~)J  (4.2.9) 

where 


V(p(u  t),x   r,r')  =  f  V(/2j3  (7-  p(u.,T)  +  x.  , 


0 


T+i  -  1 

+  T(~i^i-1})  +  I  +   (^~^)(f ''f))dT        (4.2.10) 


and 


+   1 
*i  =  2  (^3  Zi  -  *i^ 

£  (~i'~i,T)  "2  (v^3  £(yi'T)   ±  P(ui»T)>  • 

We  now  express  the  integral  over  t  as  a  sum  of  integrals  over  continuous 
functions.   Proceeding  as  before  we  fix  u. ,v.  and  partition  the  integrals 
over  t.   The  integral,  V(p(u.  ,t)  ,x.  ,  r  ,r')  ,  is  identical  to  that  for  the 
3-particle  density  matrix  (4.2.1).   The  remaining  two  integrals  require 
six  partitions: 
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T 

q 


"i,q|.    !<q<3 


L      lVi,q-3l>    3<^6 


TQ    -    0  T      -    1    . 


Now,    let  k_,...,k-   be   defined  so   that 
U  D 


k_   =  0,      k     =   7,      and 


0  <rk     <  Tk     <   '  ' '   <  Tk     <  l 
Kl  *2  K6 


The  T-interval  [0,1]  is  then  partitioned  by  the  points  i      ...t,  : 


1  6    j+1 

J  dT  V  ...  =  2        dT  V  ... 
0  j=0  T. 


Within  the  "j      "    interval,    (t,     ,Tu        ],    it   can  be   shown  that 

.        J  J  +  1 


/2B      -  - 

7— ^  p   (u.,v.,t)  =  a.   -to 


where 


+ 


+ 


(2) 


+  S^ 


afi   =   °5 "°*n    =/p/(2n)      Z      (/3    S'    '6(q+3,kJ   +  S^'fiCq,^)) 
o  jq  £=\  °  " 


q  ■   1,2,3;    j   =   1, ... ,6 


(4.2.11) 


and 


(1) 


sigti(u     ), 


S<2)   =  slgn(v.q)    . 


Thus, 


+  +     +       +  6  j+1  +  + 

;(p"(u,v,T),x',Q",Qf")    =     Z       J  V(|u~  +  U~   T|dT 

j=0    T, 


(4.2.12a) 
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where 

+    +    +      +  +    + 

v~,   =  o~  +  *~       +  Q~  +   (^-)(Q,_  -  Q") 
~  1    ~j    ~i-  1    ~       n   ~      ~ 

+      +     +        +    + 

U"  =  -  a!  +  ^~  +  ~   (Q'~  -  Q")  (4.2.12b) 

~Z       ~0     ~1    n   —       — 

and  +       +    +    +    + 

l»>!  s  /2B  xT ,  Aw?  =  uu7  -  uT    . 
~i     K  ~i   ~i   ~i   ~i-l 


It  is  to  be  noted  that  equations  (4.2.4)  have  the  same  form  as  equations 
(4.2.12).   In  both  cases,  we  must  calculate  integrals  of  the  form: 


k  k 


1  -  J  v(iHi  +  H2  Ti)dT  =  I  v(x(T>)dT  » 


V   = 


T  .  T. 

J 


X(t)    =     /u*  +  2U-U2   t  +  \]22   t2 


or 

"k  1  1 

i.        X(T)12  X(T)6 

in  the  case  of  a  Lennard- Jones  potential.   Although  this  integral  can  be 
evaluated  exactly,  the  resulting  exact  expression  is  a  rather  involved 
function  of  at  least  three  parameters: 

— y—  ,  -r     ,  and  T.  . 

Since,  in  an  actual  computation,  the  values  of  the  parameters  involved  in 
V1  would  correspond  to  a  large  number  of  Monte  Carlo  events,  it  would  be 
time  consuming  to  evaluate  the  exact  expression  each  time  V  is  needed; 
furthermore,  a  table  of  integrals  would  be  prohibitively  large.   By 
approximating  V1  by  a  quadrature  formula,  the  time  of  calculation  of  V 
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can  be  kept  small  without  sacrificing  accuracy.   The  two  point  Gauss- 

Legendre  quadrature  has  the  advantages  of  being  simple  and  having  an  error 

term  proportional  to  (t,-t.)  V  "    (as  in  Simpson's  rule).   Using  this 

k  J 

approximation  we  have, 

V'  =  A[V(t")  +  V(t+)]  +  0(A5) 


where 


T.-T.       + 
A  =   2  -1  and  t"  =  t.  +  A(l  +  .57735027) 


4.3   Strategy  for  the  Monte  Carlo  Approach 

We  now  discuss  the  general  procedure  to  be  followed  in  the  evaluation 
of  the  multiple  integrals  used  in  the  calculation  of  the  pair  distribution 
function  by  the  Monte  Carlo  method.   Equations  (2.4.4-7),  (4.2.1)  and 
(4.2.9)  suggest  that  in  general  we  will  be  concerned  with  multile  integrals 
of  the  form, 


where  the  integration  "over  the  third  particle,"  r„ ,  in  equations  (2.4.4-7) 
has  been  transformed  to  the  relative  coordinate,  Q  =  r„  -  t  (r1  +  r») . 
The  dimensionality  of  the  Q-space  integration  can  be  reduced  by  considering 
the  following  coordinate  system:   Let  us  take  r  =  r„  -  r1  along  the  Q  -axis 
as  in  Figure  2  with  the  center  of  mass  of  particles  1  and  2  at  the  origin. 
Then  Q  is  the  radius  vector.   Since  the  integrand  of  7?  is  actually  a 
function  of  two-  and  three-particle  density  matrix  elements  of  particles 
1,  2,  and  3,  and  since  the  interactions  are  spherically  symmetric,  there  is 
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azimuthal  symmetry  with  respect  to  Q,  i.e.  the  density  matrix  elements 
shall  be  independent  of  0.   Hence  in  (4.3.1)  we  are  free  to  set  0  =  — 


Figure  2.   Coordinate  System  for  Calculations 


;o  that  Q  =  (0,0  ,Q  ).   The  integral  over  Q  is  then  restricted  to  the 


V 


■Q     right  half-plane   and   takes    the    form: 


,   3n        ,   3n 
du        r   dv 


*  -   2TT  J   V%  I   dQZ  I   %n  I   ^xn  J  3T  J  T31T  NS'M*£> 


(4.3.2) 


To  convert  the  multiple  integral,  7?,  into  a  form  useful  for  Monte  Carlo 

sampling  we  must  first  consider  a  probability  distribution  for  0  ,Q  .   We 

2 
assume  for  now  that  a  suitable  distribution  function,  F(Q  ,Q  ),  has  been 

y  z 

chosen.   The  integral  then  becomes 

i  ■  1 1  «£*».  ra£v  j  %n  j  d,xn  i  if  j  ^ 

0-oo  2      2 


x  [tt 


F(Q,x,y,u,v) 


ttF, 


r(Qy,Qz) 


■]  =  <F->  • 


(4.3.3) 
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The  multiple  integral,  71,    could  at  this  point  be  estimated  by  a  Monte  Carlo 
procedure.   The  variance  can  still  be  reduced  somewhat,  however,  by  further 
importance  sampling.   It  will  turn  out  that  part  of  the  function,  F,  can  be 
factored  out,  normalized,  and  used  as  a  probability  distribution  for  the 
selection  of  certain  variables.   Let  these  variables  be  symbolized  by  the 
vector  §  and  the  rest  by  T\.      Thus  F(§,T])  =  U(5)G(|,T)),  (U(?)  >  0)  and 

71   =  J  d|w(|)  J  dT]aj(Tj)F(|,Tj)  (4.3.5) 

=  J  d§w(|)U(|)  J  dTp(T))G(|,Tp  (4.3.6) 

where  w(§)  and  ^(§)  are  weighting  functions. 
Now  define  the  function 

U(§) 

P(|)  =  -f-  (4.3.7) 

where 

Z  =  J  d|w(|)U(|)  (4.3.8) 

so  that 

7?  =  Z  J  d|w(|)P(|)  J  dTjiu(Tj)G(|,7j)  .  (4.3.9) 

The  variables  ^  a*e    thus  weighted  by  the  probability  distribution  w(§)P(§). 
It  is  intuitively  apparent  that  it  is  advantageous  to  use  P(§)  in  the 
probability  distribution  for  ^  since  it  biases  samples  of  §  in  favor  of 
those  for  which  the  factor  U(^)  in  F(^  ,Tj)  is  relatively  large. 

In  general,  it  is  difficult  if  not  impossible  to  generate  random 
variables  with  distribution  P(?)w(^)  directly.   Instead  of  taking  independent 
samples,  it  will  prove  easier  to  generate  a  Markov  chain  of  random  variables 
£(k)  whose  distribution  converges  to  P(§)w(^)  as  k  -»  co .   In  the  case  of  a 
discrete  sample  space,  Metropolis  et  al.  [30]  showed  how  to  do  this  in 
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calculating  averages  with  the  classical  Boltzmann  distribution,  and  the 
method  was  formalized  by  Hammersley  [31].   Jordan  [9]  applied  the  method 
to  the  evaluation  of  Wiener  integrals  by  extending  Hammersley' s  argument  to 
the  case  of  a  continuous  sample  space.   The  present  approach  stands  some- 
where between  that  of  Hammersley  and  that  of  Jordan. 

Let  us  assume  that  in  a  Monte  Carlo  scheme,  the  first  choice 
§ (0)  =  ^  is  determined  by  some  arbitrary  distribution  (e.g.  if  it  is 
always  the  same  number,  the  distribution  of  |(0)  would  be  a  6  function). 
The  next  choice  £,(1)    =  5'  would  be  determined  by  a  Markov  transitional 
probability  density  p(§'|5)  which  satisfies  the  conditions 

J  P(|'||>d5'  =  l  (4.3.10) 

and 

w(|')P(|')  =  J  p(|'|?)p(pw(pd|   •  (4.3.11) 

Equation  (4.3.10)  is  the  usual  normalization  condition  for  conditional 
probabilities.   Equation  (4.3.11)  prescribes  that  if  §  has  distribution 
w(^)P(£)>  then  the  new  random  vector  ^',  determined  by  the  Markov  transition 
probability  p(§'|§),  will  have  the  same  distribution.   The  next  random 
vector  %(2)    ■  5"  would  be  determined  by  the  same  transitional  distribution, 
p(§"|C')  and  so  on.   The  fact  that  the  distribution  for  §"  depends  on  §' 
and  not  on  ^  is  the  Markov  property  of  the  process.   Continuing  in  this  way, 
we  would  obtain  a  sequence  of  vectors,  ^(k),  the  distributions  of  which 
converge  to  P(?)w(^),  equations  (4.3.10)  and  (4.3.11)  being  necessary  and 
sufficient  conditions  for  this  convergence. 


Another  condition,  that  the  states  §  be  ergodic  is  also  required  and  holds 
in  the  cases  studied  here. 
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Furthermore,  if  equations  (4.3.10)  and  (4.3.11)  hold  and  the  states 

%(k)   and  T|(k)  form  an  ergodic  set  and  if  71  =   J  d§w(§)P(^)  !'  dTp(Tl)G(5 ,11) 

t  M 
is  estimated  by  —  2  G[^  (k)  ,T|(k)]  where  each  Tj (k)  is  distributed  in- 
M  k=1   ~ 

dependently  by  W(\\)    (^(T])  can  also  be  regarded  as  a  Markov  transition 
probability  density  for  T|  since  it  satisfies  conditions  similar  to 
(4.3.10-11)),  then  by  the  law  of  large  numbers  for  Markov  chains  [32]  the 
estimate  will  converge  to  the  integral,  71  >   with  probability  one  as  M  -•  «, 
In  order  to  find  a  transition  probability  satisfying  (4.3.10-11)  we 
introduce  the  "trial  transition  density" 

P*(|'||)  =  w(|')p**(|',|)  (4.3.12) 

where  p**(?'>?)  is  a  function  satisfying 

P**(|',p  =  P**(|»|')  (4.3.13) 

and 

J  w<5,)P**(§,f|)d|'  =  1  •  (4.3.14) 

We  then  have  the  following 

■k 

Theorem: 
The  choice, 

HV)  p(5') 

p*<i'li>P(fr  Lfwr<1 

P*(|'||)  +  6  <|'-|)       J  d§" 
P(|'|?)=  HZ")  (4.3.15) 

^Pcfr^ 

P(S")  P(S') 

X  P*(|"||)(l  -  p^-)      if  p(f)->  l 


/V 

This  theorem  resembles  that  of  Jordan  [9].   The  choice  of  p*(§'|§)  and 
the  proof,  however,  follow  more  closely  the  approach  of  Hammersley  and 
Handscomb  [31]. 


45 


satisfies  conditions  (4.3.10)  and  (4.3.11). 
Proof: 

(1)  We  show  I  =  J  d§'p(§'|§)  =  1:   Let 

P<§') 

r(i,)H7(fr; 

then  from  (4.3.15)  and  (4.3.12) 

i  =  J  "OpP^CI'^Mpdi'  +J  w(|')p**(|'|)d5' 
{|':r(|')  <  1}  U':r?')  >   l) 

+  J  w(5")p**(|",|)(l-r(|"))d|" 
U":r(|")  <  1} 

the  first  and  the  second  half  of  the  third  terms  cancel  leaving 

I  =  j  w(|,)p**(|,,|)d|'  =  1  . 

(2)  We  prove  (4.3.11)  first  by  showing 

p(|'||)P(|)w(|)  =  P(?||')p(i')w0p  •  (4.3.16) 

Case  1:   If  P(§')  <  P(£): 

PCS1) 

p(|'U)p(|)w(|)  =  p,v*(r»|)w(|*)  p^y-  p(|)w(|) 

=  P**(|,|')w(|)P(|')w(|') 
-  p(§|?,)P(?')w(?')  • 
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Case   2:      If  P(§')  >  P(§),  §'    *  §:    (if  ?'    =  ?    (4.3.16)    holds   trivially) 

p(|'||)P(|)w(|)   =   p**(|',|)w(|')P(|)w(|) 
=  P**(|,|*)w(|)P(|)w(|') 

PCS) 
=  p*(|»|')  pfry  p(|')w(|') 

=  pCglpPCpwCi')  . 

Thus  (4.3.16)  holds  in  all  cases.   Hence 

J  p(ripp(|)w(pdi =  j  p(iir)p<r)w(s,)di 

=  P(|*)w(|')  J  P(|||')d|  =  P(|')w(|') 

Q.E.D. 

In  the  application  of  the  above  theory  to  the  calculation  of  the  pair 
distribution  function  terms  n^   ,  n   ,  n_  ;  (equations  (2.4.3-7)  the  vector 
§  will  have  components  which  are  path  break  points  (e.g.  x) ,  "Simpson's 

rule"  parameters  (e.g.  u)  and  (in  the  case  of  n^;   (r)  and  nj   (r))  the 

2 
coordinates  of  the  third  particle  0  ,  Q  .   In  these  calculations, 

P**(§>?)  is  chosen  to  be  a  function  of  the  path  variables  only.   Further, 

in  the  case  that  ^  involves  two  paths  x  and  y  (as  with  np   (r))  we  assume 

~  2 

the  "trial"  paths  to  be  chosen  independently  with  the  same  distribution, 

i.e. 

P**(J',|)  =  f(x',x)f(y',y) 

where  x  and  y  are  3n-tuples  representing  the  path  break  points.   The 
"trial  transition  probability"  for  the  joint  distribution  then  factors: 
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p*(|'|p  =  f^,^)wi(^,)f(y,»y)wi(y)--- 

where  w. (x)  is  the  weight  function  for  the  path  x  and  is  given  by  the 
relation 

w.(x)dx3(n"1)  s  du.   .  (4.3.17) 

1  ~  rnx 

We  determine  f(x',x)  by  requiring  for  simplicity  that  the  "trial 
transition  probability"  be  Gaussian.   We  thus  choose  the  form  for  f(x',x) 

n  2 

w1(x)f(x',x)  =  C  exp[-  |  Z      (ax^  +  px^  +  yx     +  &x   L>  ] 

i=l 

3n 

=  (2tt)3/2(^)2  exp[-|  Z   (x1-xi_1)2]f(x',x)     (4.3.18) 

i=l 

where  we  have  used  equations  (4.2.2)  and  (3.1.10)  in  expressing  w..(x). 
The  symmetry  requirement  (4.3.13)  implies  f (x' ,x)  =  f(x,x')  which  from 
(4.3.18)  yields  the  conditions: 

2  .     2 
a  -  1  =  y 

2        2 
P  -  1  -  6 

ap  +  1  =  Y6 
a6     -  PY  • 

The  second  of  these  equations  is  implied  by  the  first,  third,  and  fourth 
equations  so  that  only  three  are  independent,  and  a,  (3,  and  Y  can  be 
determined  in  terms  of  the  arbitrary  constant  6.   The  normalization 
condition  (4.3.14)  determines  C  in  terms  of  this  same  constant.   We  thus 
finally  obtain: 
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3n 
2,    2 


/    xir/    i      \        /o   \3/2    ,n(l+6   K  /  ,,  ,.2.    n 

w1(x)f(x',x)   =    (2rr)  ( — ^ — «■)        exp{-(l+6   )  -^ 

X     S     [x!-x.    .    --~=7   (x  -x     J]2).  f4.3.19) 

i=l     -1  -1"1        /i+6       ~     ~ 

If  we  make    the    transformation  to   the   3n-tuple,    x  =   (x, ,x2,...,x  ): 

x'    =  —=-   (X  +  6x) 
~        /ST2        ~  ~ 

then  X  has  the  distribution  w  (x)  and  x  =  X-  =  0- 

In  summary,  the  random  sequence  (^  (k)  ,T|(k))  is  generated  by  the 
following  iterative  procedure.  Assuming  we  have  chosen  %,  (k-1)  ,T^(k-l)  and 
have  determined  U(§(k-1)): 

(1)  Choose  Tl(k)  according  to  the  distribution  w(T|)  . 

(2)  Choose  the  "trial"  vector  §'  in  the  following  way: 

Paths:   The  3-vectors  x-ijXoj'-'jX   -i  are  generated  by  the  weight 
function,  w. (x) ,  i.e.  according  to  the  "interpolation"  formula 
(3.1.9)  (x  =  X  =0)-   The  "trial"  path  coordinate  is  then 
determined  by  the  formula: 

x!  =   *    (x,  +  6x.(k-l))     i  =  0,1,. ..,n  .    (4.3.20) 
~x    /l+6 

Other  variables  in  £ ' :   "Trial"  choices  for  variables  which  are 


not  path  variables  are  made  according  to  their  corresponding 

weight  functions. 

P(S')       U(§') 
(3)   The  ratio,  r(|')  -  p(^(k-l))  =  u(%(k-l))  '  is  calculated-   If 

r(§')  >  1  then  5(k)  -  §*.   If  r(§')  <  1  then 
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r    I'  with  probability  r(^') 
|(k)  =/ 

^  5(k-l)  with  probability  l-r(§')  . 

The  problem  of  choosing  the  initial  vector  £,(0)    still  remains.   In 
theory,  any  choice  will  yield  a  sequence  converging  to  a  random  variable 
which  has  distribution,  P(^)w(§).   The  number  of  steps  required  for  con- 
vergence depends  (in  the  probabilistic  sense)  on  the  distribution  of  the 
initial  vector,  ^(0)-   The  selection  process  will  be  discussed  later, 
however . 

n 

The  multiple  integral,  —,    is  thus  estimated  by  the  sum, 

1   M 

—  E  G(^  (k)  ,T](k)) .   The  normalization  function,  Z,  can  be  estimated  by 

M  k=l   ~ 

Monte  Carlo  sampling.   If  6  =  0,  the  estimate  for  Z  comes  out  as  a  by- 
product: 

1  M 
Z  =  ±     Z     UO-'Ck))  (4.3.21) 

k=l  ~ 

where  Z  '  (k)  are  the  "trial"  vectors  used  to  determine  §(k). 

Another  advantage  of  taking  6=0  arises  when  the  same  random  numbers, 
and,  therefore,  the  same  trial  vectors,  are  used  for  calculations  involving 
only  different  values  of  r.   Consider  a  particular  event,  k,  in  a  Monte 
Carlo  sampling  process.   Assume  a  trial  vector,  5'>  has  been  selected  and 
is  tested  as  in  step  (3)  above.   Whether  it  "succeeds"  or  not  depends  on 
its  value  and  the  value  of  r.   Suppose  that  the  calculations  are  performed 
for  a  set  of  values  of  r.   There  will  in  general  be  a  group  of  calculations 
over  a  set  of  r  values  such  that  the  trial  vector  succeeds  and  ^(k)  =  ?'• 
There  will  be  other  groups  involving  r  values  such  that  the  appropriate 
vectors  are  given  by  perhaps  ^(k-1),  or  perhaps  §(k-2)  and  so  on.   Each 
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group  involving  a  different  value  for  % (k)  will  in  general  refer  to  a  sub- 
set of  one  or  more  r  values.   By  carefully  keeping  track  of  these  groups 
and  performing  calculations  which  are  independent  of  r  only  once  for  each 
group,  considerable  computational  time  can  be  saved. 
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5.      EVALUATION  OF   THE   PAIR   DISTRIBUTION   FUNCTION   TERMS 

5.1      The   Cyclic   Exchange   Term 

The   simplest    three-particle   contribution  to   the   pair   distribution 

function    (2.2.6)    is    the   cyclic   term   (2.4.5).     As   suggested   in   the   previous 

section,   by  applying  equations    (4.1.8)    and    (4.2.5)    to  equation    (2.4.5)   we 

obtain 

3        2  Q* 

n£1}(r)    =   2e      4p  J   d3Qe  P     E{f[T)(t)  ,r  ,Q_] 

X   F[l]"(T),Q",-Q+]F[T]+(T),Q+,-r]|T)(l)    =    q(l)    =    0} 

(5.1.1) 

+ 
where   F,   T]    (t)   were   defined   in    (4.2.6-8).      The  multiple    integral  approxi- 
mation  is    then  employed,    and,    as    in  equation    (4.3.1),   we   have 


yn  -    2' 


where 


*  l$r  «*- 5  .=,  I  (v"  +  v23 +  »S>*1 +  "(-S) 

2  i=l   o  n 


(5.1.2) 


v"  =  V(/23(t-  P(u't>   +  x-,-    i   +  t(x.-x,    ,))■  +  r  + 

iz  yn  «,   ~  ~i-i  ~i   ~i-i  ~ 


(I±rli)(Q"-r))  (5.1.3a) 


V3J  2  V(/23(t-  P"(u,v,t)   +  x"        +  T(x"-x7    .))  +  Q" 
zj  yn  ~      ~   ~  ~i-i  ~i   ~i-i  ~ 


-    2    (I±n^i)    Q)  (5.1.3b) 
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V2*  h  V(/2p(7-  p+(u,v,T)  +  x|   +  T(x+-x+  ,))  +  Q+ 


(I:T:i)(S+  +  5)}  '         (5.1.3c) 


Consider  now,  the  integration  over  the  third  particle.   According  to  the 
discussion  in  Section  A. 3: 

Jd3Q  ...  =  tt  J  dCyiQz  ...  where  Q  =  (0,  Qy,  Qz>  . 

Choosing  spherical  polar  coordinates,  Q  =  Q|sin  9|,  Q  =  Q  cos  6  and 

Q  =  |q|  We  get 

1 
J  d  Q  ...  =  2tt  J  Q  dQ  j  d(cos  6)  ... 
0     -1 

2 
2  -0  /B  2  2 

The  factor  Q  e    K  suggests  that  we  choose  0  according  to  a  ^ 

distribution  with  three  degrees  of  freedom: 

-I  -1     1    2 

T(Q2)  =  2tt  2  B   2  (Q2)2  e"Q  /P  (5.1.4) 

J  T(Q2)dQ2  =  1  . 
Equation  (5.1.2)  then  becomes: 


n(»(t).£e"*er2jr«wJ*<£pl 

/2         0         -1 


_____  Q  O'i 

where  X  =  /2nB  is  the  thermal  wavelength.   The  factor,  exp[-  *-  Z  V   dT] 

n  i = 1  «   12 
1  L  0 

can  be  taken  as  a  weighting  function  for  importance  sampling  with  a  Markov 
process.   We  therefore  define  the  probability  distribution: 
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n   1 
Pc(x,u,Q2,cos  6)  =  [Zc(r)]_1  exp[-  |  S  J  v"  dT]     (5.1.5) 

i=1  0 


where  the  normalization  condition  requires: 


0-1  2 


X  exp[-  |  2  J  V^2  dT]  .  (5.1.6) 

1=1  o 


The  multiple  integral  finally  takes  the  form: 

3   -^-r2 
nC1>(r)  =  7"  e  4?    ZC(r)8C(r)  (5.1.7) 

where 

3n         n  l 
«C(r)  2  J  dPC  S   %n  I  3T  «*rf-  I  *  I  ^Xs^   (5-1-8) 

and  we  have  used  the  shorthand  notation: 

J  dPc  ...  »  J  r«W  J  f  J  *   J  if  Pc(x,u,Q2,a)... 

0-1  2 

dv3 
where  a  =  cos  9.   Since  |  dP   [*  d|jl   P  — ^—  =  1,  gr(r)  can  be  regarded  as 

a  weighted  average  and  can  be  estimated  by  Monte  Carlo  sampling  as  outlined 

in  Section  4.3.   Thus 

M  n   1 

gc(r)  -  ^  I  exp[-^  2  J  (V^(k)  +  V^(k))dT]       (5.1.9) 
k=l        1=1  0 
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where  the  argument  (k)  implies  that  the  variables  of  integration  are  now 

2 
the  realization  of  a  stochastic  process:  x(k) ,  y(k),  u(k),  v(k),  Q  Ck), 

2 
a(k)  where  x(k),  u(k) ,  Q  (k) ,  and  a(k)  are  components  of  a  Markov  process  and 

y(k)  and  v(k)  have  independent  steps. 

5.2   The  Direct  Term 

5.2.1   Separation  of  the  "Hard-Core"  and  "Correlation"  Terms 
The  three-particle  contribution  to  the  pair  distribution  function 
which  does  not  depend  on  the  exchange  of  particles  will  be  referred  to  as 
the  "direct  term,"  n   (r) ,  and  is  defined  by  equation  (2.2.4).   Choosing 
our  usual  center  of  mass  coordinate  system  and  making  use  of  the  azimuthal 
symmetry,  the  integral  becomes 

CO         CO 

n^  =   tt  J  dQ2  J  dQz  [W3-W(W+  +  W"-1)]  (5.2.1) 

where  we  use  the  shorthand  notation: 


W3  =  \y<l  2  3|e   3|l  2  3)  =  (e    12   2J   13  > 

r  -pH  -p  V 

W  =  X6(l  2|e   2|l  2)  =  <e   12> 


W+=  X6<1  3|e   2|1  3)  =  <e    13> 

,       -PH  -PV 

W'  =  Xb<2  3|e   Z|2  3>  =  <e   2J) 


(5.2.2) 


The  expressions  on  the  extreme  right  side  of  (5.2.2)  represent  path  averages 
and  V. .  is  the  "time  averaged"  interaction  between  two  particles  whose 
initial  locations  are  r.  and  r.  respectively  (see  Section  3.2). 
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It  is  seen  by  equation  (5.2.1)  that  one  difficulty  in  the  calculation 

of  the  direct  term  is  that  there  is  no  obvious  weighting  function  for  the 

»   3 
distribution  in  Q  space.   The  integrand  of   d  Q  ...  goes   to  zero  for  large 

|q|  by  virtue  of  a  cancellation  between  the  two-  and  three-particle  density 

matrix  elements.   As  a  first  step  in  considering  the  integration  over  Q, 

we  make  use  of  the  apparent  symmetry  of  the  integrand  with  respect  to  the 

x-y  plane: 

CO         CO  CO         CO 

tt  J  dQ~  J  dQ  [W  -W(W  +W"-1)]  =  2rr  J  dQ  J  dQ  [...] 
0    -co  0   y  0 

so  that  integration  need  only  be  considered  in  the  first  quadrant.   Let  us 
call  this  quadrant,  [0  >  0,  Q  >  0],  S.   Then  let  us  partition  S  by  a 
circle  of  radius  Rn  about  particle  2.   The  set  S  is  then  considered  to  be 
made  up  of  two  subsets  Sn  and  S'  where 

so  =  *V  V  %  >0>   QZ  >  °'  Qy  +  (QZ  "  2  r)2  *  %] 

and  S1  =  compliment  of  S   in  S  =  S/Sn  =  S^.. 

Thus  we  may  express  the  integral  over  Q  as 

CO         CO 

2TTJ  di!  d%  •••  =  2ttJT  did%  ••• 

o         o  s 

-  2*n  did*z  ••• +  2n n  dQy  dQz  •••     (5-2-3) 

sQ  s> 

Recalling  that  the  interaction  potential,  V(r),  is  positive  for  r  <  1  and 
increasing  rapidly  as  r  becomes  smaller  we  expect  that  the  average  inter- 
action, V..,  should  be  quite  large  for  small  values  of  r.-r.  <  1.  Thus, 
we  expect  that  there  should  be  some  maximum  value  of  L,  0  <  R„  <  1,  such 
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-0V23 
that  e      is  negligible  for  Q  over  the  region  S   for  most  paths  so  that 

-0V12   -pV23   -pv13        -£V23 
(e      e      e     )  and  (e     )  can  be  neglected.   We  will  call  Rn 

the  hard  core  radius  and  S   the  hard  core  region  since  it  is  in  this 

region  that  the  hard  core  (strongly  repulsive)  interaction  between  particles 

2  and  3  dominates  the  density  matrix  elements.   Having  chosen  such  an  R,  , 

the  direct  term  becomes: 


n£1}(r)  =  2tt  W  JJ  dQ2  dQz(l-W+) 


y 
so 


+  2tt  JJ  dQ  dQ  [W  -W(W++W"-1)]  (5.2.4a) 

S* 

=  n_  (r)  +  n_,(r)  .  (5.2.4b) 

The  first  integral  involves  only  a  two  particle  term,  W  ,  and  can  be 
evaluated  by  quadrature  over  a  table  of  values  for  W  .   The  details  will 
be  described  in  Appendix  A.   The  second  term  can  be  expressed  as  one  Wiener 
integral: 


2  Jn   ,  -^12,  "PV23  -PV13    .    -PV23,    ,  "PV13. 
I      dQ 
y   z 


2tt  JJ  dqj  dQz  (e   12(e   23e    13  -  <e   23>  -  (e    13>  +  1)) 


-pv23  -pv23 

where  in  the  second  term  the  (e     )  cannot  be  replaced  by  e      since 

-pv12     -pv23 

there  is  a  correlation  between  e      and  e      so  that 

-pv   -pv      -pv    -pv  -pv   -pv 

(e        ll  e        23>   f   <e        12><e        23>,    and   similarly,    <e        U  e        13>   ^ 

-pv12   -pv13 

(e     )(e     ).   It  will  prove  useful,  however,  to  employ  this 
approximation  and  split  the  integral  in  the  following  way: 
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1  2 

ngl  (r)    -     ns,  (r)   +     tig,  (r) 


whe  re 


and 


i  ..       2  "pvi2     "pv23  _Pvn 

ns,(r)    =   2tt  JJ   dQ^  dQz<e        L/(e  -l)(e        iJ-l))  (5.2.5a) 

S' 


2     r  .     2    rr  ho2  ho  i  "pVl2r  ~pV23  <  "pv^ 

ns,(r)  =  2tt  JJ  dQ  dQz<e     (e     -<e     ) 


y 
s' 


-pv     -pv 

+  e    1J  -  (e    13))>  .      (5.2.5b) 


That  the  two  integrals  should  converge  is  shown  by  the  following  argument. 
Let  us  consider  the  asymptotic  behavior  of  the  integrand  of  nQ,(r). 
According  to  our  interpretation  in  Section  3.2,  V..  is  the  time  averaged 
interaction  between  particles  i  and  j  whose  uncertainty  in  position  is  of 
order  \  =/2n(3.  As  |Q|  -*  «  we  expect  V   and  V   to  approach  zero  like 

1  5 

9 (—7)  so  long  as  JO  J  »  CTMAy,  the  maximum  deviation  of  the  particles  from 
locations  1,  2,  and  3  for  a  particular  path.   According  to  the  conditional 
Brownian  measure,  these  paths  which  deviate  from  their  initial  positions 
a  distance  more  than  X  do  so  with  exponentially  small  probability.   Thus, 
for  each  X  there  is  a  maximum  distance,  IL.  ,  such  that  if  |Q  |  >  RwAy 
then  |Q  j  »  (TwAy  for  "most"  paths  and  V   ~  0(— 7-),  i  =  1,2.   Now  under 


r13 

these  conditions 


-pv        -PV  _ 

(e      -  l)(e      -  1)  =  P  V23V13  + 


13  23 


for  "most"  paths.   Thus  we  expect 


<e    12(e  '  23  -  l)(e    13  -  1))  -  O(-i-) 

Q 


as  Q  -»  co.  The  integral,  n  ,(r),  therefore,  is  expected  to  exist.  Further- 
more the  above  argument  applies  if  the  integrand  is  replaced  by  its  absolute 
magnitude  so  that  by  Fubini's  theorem,  the  integral  with  respect  to  the 

product  measure  of  Wiener  measure  and  ordinary  Lebesgue  (Q)  measure  exists. 

2 
A  similar  argument  applies  to  n  , (r) . 

2 
As  was  mentioned  before,   n„,(r)  depends  upon  the  correlation  between 

-pv12         -ev23        -pv13 

e     ,  and  e      and  e     .   In  fact 


ns,(r)  =  2rr  W  JJ  dq^   dQz[W  p(e 

S'  • 

-   -^V12   "PV23 
+  W  p(e     ,  e   ZJ)]  (5.2.6) 

where  p(x,y)  =  -^ )  (/(  ^ '  ' '    is  the  correlation  coefficient  of  x  and  y 

with  respect  to  conditional  Wiener  measure.  As  r  and/or  Q  becomes  large, 

-pv12  -pv23      -pv13 

e      and/or  e   _  ,  e      approach  the  constant  unity  for  "most"  paths. 

-pv12 

Also  as  r  -»  0,  e      -*  0  due  to  the  repulsive  core,  and,  in  either  case, 

2  2 

n  ,  (r)  -»  0.   Hence  we  expect  that  the  correlation  term  n  (r)  need  be 

calculated  for  only  a  few  values  of  r.   Let  us  now  express  the  terms 

1  2 

n.,(r)  and  n„,(r)  in  the  multiple  integral  approximation  using  the  center 

of  mass  transformation  and  equations  (4.2.1)  and  (4.2.5).   We  have 

3n 


Lns,(r)  =  2tt  JJ  <  dQz  J  d,xn  J  d,   J 


du" 


y   z  J   rxn  J   ryn  J  93n 
S'  l 


,  3n 
X  I  *~-   D(D  -D(D'-l)  +  0(^) 
2  n 


(5.2.7) 
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A     3tl 

du 


-ns,(r)    -   2n  JJ   •%   dQz  J  d,^  J   d%  J  -fj 


whe  re 


8<  2" 

,   3n 
X  J  ~p  D(D   -W  +D"-W")    +   O(^)  (5.2.8) 

2  n 


n      1 


S.   v    r  \^/"?q  c-L 


D  =   D(x,u,r)    =  exp[-  ?■     Z       '   V(/2p    (-—  p(u.  ,t)   +   x. 

~    ~  n    .     ,    «J  Vn    ~    ~1  -1-1 


1=1 


0 


/n 
+  T(xi-x._1))    +  r)dT]  (5.2.9) 

+  +  +  a       n       l  _       i  + 

D     h   d    (x,y,u,v,Q   )    =  exp[-  *     Z       '   V(/2p    (—   p    (u,  ,v.  ,t) 
n   i=l  J0  /n  ~      ~L   -1 

+  +  +  + 

+  xT    ,    +  t(xT    -    xT    ,))   +  Q")dT]    .  (5.2.10) 

~i-l  ~i        ~i-l  „  J 


The  factor,  D,  appears  in  both  integrals  and  can  be  used  in  importance 
sampling  with  respect  to  u,x.   Hence  we  define 

/_n  D(x,u,r) 

w 
where  the  normalization  factor 

is  an  approximation  (to  order  —z)    to  the  two-particle  Boltzmann  Slater  sum, 

n 
W.   The  integrals  (5.2.7-8)  then  become 


\,M    -  2n  W(")(r)  It   dQ2dQz  J  dP<n)  J  d,, 


y   z  o    D   J    yn 
S' 

,  3n  , 

X  f  ^—  (D  -1)(D~-1)  +  O(^)    (5.2.12) 
J  23n  n2 


S' 


X  (D+-W++d"-W~)  +  0(-^)      (5.2.13) 


where  we  have  used  the  notation, 

dp(n)  a  p(n)      }  du__ 

D      D   ~'~'    23n   *xn 

These  integrals  can  be  evaluated  by  Monte  Carlo  sampling  in  the  product 
space  of  all  the  variables,  x,y,u,v,Q  ,Q  ,  provided  an  appropriate  prob- 

~  ~  ~  ~  y  z 

ability  distribution  for  Q  is  selected. 

To  get  an  idea  of  the  behavior  of  the  integrands  as  a  function  of 
Q,  the  random  function,  F*  =  D  (x(k) ,y (k) ,u(k) ,v(k) ,Q  ),  has  been  calculated 
by  controlling  the  Q  values  and  letting  the  other  variables  vary  randomly 
according  to  their  various  distributions  for  each  Q.   In  Figure  3  and  4  we  plot 
for   two    temperatures,  F*  for  values  of  Q  along  a  diagonal  line  of 
length  R  which  ends  at  the  0  -axis.   The  apparent  randomness  is  due  to  the 
variation  of  the  "path"  variables  u,v,x,y.   We  see  that  for  higher  temper- 
atures there  is  a  definite  overall  structure  so  that  the  Q  distribution  is 
important.   But  as  the  temperature  goes  down,  the  variation  of  F*  due  to 
the  "path"  variable  is  of  the  order  of  the  overall  change  with  respect  to 
Q  ,  and  the  structure  of  F*  with  respect  to  Q  becomes  less  obvious,  there 
being  only  a  tendency  towards  zero  for  large  | Q~ | . 

Thus  in  considering  a  distribution  for  Q  for  low  temperatures  it  is 
not  necessary  or  even  advisable  to  devise  a  precise  function  of  Q.   Instead, 
it  will  be  sufficient  to  partition  the  Q-space  in  accordance  with  the  general 
tendency  of  F*.   Further,  we  will  assume  that  the  structure  of  F*  depends 
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2.80  3.20  3.60 


4.00 


Figure   3.      Plots    of   F*  vs.    I Q~ I     for   T  =    5°K,    r   =    1 


62 


|Q"I 


<wryw»^^pww^»»'^i»» i^"i< 


rsr 


3.20 


3.60  q.oo 


Figure  4.      Plots   of  F*  vs.      Q        for  T=5K,    r=1.9 
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Figure  5.   Plot  of  F*  vs.  |Q  |  for  T  =  1.7  K, 
r  =  1.9. 
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mainly  on  the  absolute  distance  between  particles  2  and  3,  Q  .  Any 
influence  of  Q  would  be  due  to  the  correlation,  via  the  path  T)(t),  with 
particle  1.   The  validity  of  this  assumption  can  be  demonstrated  by 
comparing  F*  for  r  =  1  and  r  =  1.9  in  Figures  3  and  4.   The  same  general 
dependence  on  Q  is  observed  for  the  different  r  values.   Furthermore  the 

exact  same  argument  applies  to  D  which  involves  the  interaction  between 

+  + 

particles  1  and  3  so  that  D  can  be  regarded  as  strongly  influenced  by  Q  . 

5.2.2  Stratified  Q-Space  Sampling 

We  now  turn  our  attention  to  the  term  n  , (r)  defined  by  equation 

(5.2.12).   Although  the  same  procedure  to  be  outlined  here  and  in  part  C 

2 
applies  to  the  calculation  of  n^Cr)  defined  by  (5.2.13),  it  has  not 

2 
yielded  satisfactory  results.   The  calculation  of  n  , (r)  will  be  discussed 

in  Section  5.2.4. 

The  integral,   n  ,(r),  is  of  the  form 


^g.Cr)  =  2tt  W  JJ  dQ^  dQz  J  dP(D+(|,Q+)  -  1) 


y 
s? 

x  (d"(|,q")  -  i)  +  oo^) 

n 

where  §  represents  the  variables  u,v,x,y,  and  dP  is  the  measure  of  £• 

+  + 

Since  D  depends  predominantly  on  Q  and  D  depends  predominatly  on  Q  , 

it  will  be  useful  to  transform  to  the  Q  ,Q  coordinate  system: 

-  2    2         12 

(q  )    -  q;  +  (Qz±|  r)Z  . 

The  integral  then  becomes 

Xns,(r)  =  71  W  JJ  Q+dQ+Q"dQ"  J  dP(.  ..)(...)  +  0  (^)     (5.2.14) 
oi  n 
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where  S'  is  now  a  set  in  the  Q  ,Q  quadrant,  S  (we  use  the  same  notation 
for  convenience):  S'  =  [Q  ,Q~:  Q~  >  RQ]  (IS  where 
S  =  [Q+,Q":  Q+  >  Q"  >  0;  Q+  +  Q"  >  r  >  Q+  -  Q-}. 

Let  us  approximate   n  , (r)  by  an  integral  over  the  truncated  region, 

S"  =  [Q+,Q~:  r°  <  q"  <  r„]  n  s  , 


where  it  is  assumed  that  R.  is  large  enough  so  that  the  truncation  error 
is  negligible.   We  then  partition  this  region  by  the  sets,  S.,  i  =  1,...,7 

formed  by  the  interactions  of  pieces  of  concentric  rings  of  radii, 

t  +         -  t 

R,  ,  R„,  R  •  S.  =  (Q  >  Q  :  R.  >  Q  >  R.  n 1  OS,  i  =  1,2,3  which  will  be 
1'   2'   3   i   L         i  —    —  i-lJ 

called  "basic  sets."   We  let 

s.  =  s.  n  s~         s.  =  s_  n  s7 

111  4    J     1 

s„  =  st  n  si         s.  =  st  n  s~ 


2  2     1  5    3    2 

3  "  S2  °  S2  S6  "  S3  "  S3 


s7  =  s-  u  s;  u  s;  U  S+ 

where  S   is  the  compliment  of  S  with  respect  to  S .   It  is  easy  to  verify 

7 

that  S.  OS.  =0  (i^i)  and  U   S.  =  S".   An  example  is  illustrated  in 
l    j  J      i=l   i 

Figure  6.   The  sets,  S.,  have  been  chosen  so  that  within  each  region, 
(D  (?>Q  )(D  (§,Q  )  -  1)  does  not  change  appreciably.   Thus,  in  the  example 
of  Figure  3  if  we  choose  perhaps  R_  =  .8,  R,  =  1.6,  R„  =  2.4  and  R„  =  3.2, 
there  would  be  relatively  little  overall  change  of  F*  with  respect  to  Q 
between  these  radii.   Likewise  1-D  (Q  ,^)  should  have  little  structure  with 
respect  to  Q  within  the  same  intervals.   Hence  the  product 
(D  (c,Q  )  -  1)  (D  (£,Q  )  -  1)  will  not  change  appreciably  over  the  inter- 
sections of  "basic  sets." 


Figure  6.   Example  of  Region  Partitions  for  S": 
The  Regions  S.,...,S_ 


In  summary,  the  integral  in  equation  (5.2.14)  over  the  infinite 
region,  S!,  is  approximated  by  an  integral  over  the  finite  region,  S" , 
which  is  partitioned  into  subregions,  S.,  and  we  have 


Vj.Cr)  ~   LnSI1(r)  =  2  W  E  JJ  d(Q+)2  d(Q~) 


7 

Z 

i=l 


S. 

l 


.+  ..+ 


X  J  dP(D  (Q  ,|)-1)(D"(Q",|)-1)   (5.2.15) 


where 


+ 

q" 

% 


-  (0,Qy,Qz±2  -> 


h^2 
(Q+)2 


(Q')2] 
(Qz+|r)2 


(5.2.16) 


J 
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Since  it  is  assumed  that  there  is  no  real  structure  of  the  integrand  with 

respect  to  Q  ,  Q   in  each  region,  it  will  be  reasonable  to  assume  uniform 

■K  2        -2 

distributions  of  (Q  )   and  (Q  )   over  the  regions  S.  for  Monte  Carlo 

sampling.   Thus,  we  estimate  (5.2.15)  by 

7  Mi 

1=1        1  k_i 

X  (D'O^i'lk)"1)        (5.2.17) 


+.2 


7 


where  A(S.),  the  normalizing  constant,  =  ff  d(0  )   d(Q  )  ,  and  M  =  S  M. . 
i  J  J  i=1   i 

bi 
Equation  (5.2.17)  suggests  that  the  same  random  variables  §   are  used  for 
each  region  S..   Although  this  simplification  is  justified  in  the  sense 
that  the  estimate  nQM(r,M)  -•  nQ,,(r)  as  M  -*  »  (provided  M.  -*  «,  VS.  f   0), 

correlation  between  the  regional  estimates  over  each  S.  does  slow  down  the 

°  l 

overall  convergence.   It  is  assumed,  however,  that  the  time  saved  in 

computation  makes  up  for  this  fact.   The  method  of  generating  the 

variables  §,  was  discussed  in  Section  4.3.   Details  such  as  the  generation 

+     ~k 
of  Q  random  variables  over  the  various  regions  and  the  determination  of 

A(S.)  are  covered  in  Appendix  B. 

It  remains  to  determine  the  sample  sizes,  M.  .   We  would  like  to 

choose  them  so  as  to  minimize  the  total  variance  of  the  estimation, 

1-  7 

fi  ,,(r,M),    subject   to   the   condition,      Z     M.    =  M,   where  M.    =   0   if   S.    =0. 
b  i=l      l  ii 

We   consider   the   general   problem  of  evaluating   integrals   of   the    form: 
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g  =  J      dQ  J   dpi   f  (|,Q)   where  J*   dp,  -   1 
S" 
K 

=  E  A(V  mttJ1  d?  J" d*  f(i'2> 

J-1  J     s 


J 


where 


K 
S"    = 


U     S  ';        S.    OS.    =0,    i^j;      andA(S.)=j       dQ      . 

J-1  s. 


SJ 


Then  let  the  total  Monte  Carlo  estimate  be  given  by 


K 
L  =  2  A(S  )  <f>  (5.2.18) 


where 


<f>M.  "I"  j?  "Sk-i^' 

J    J  k=l 


(5.2.19) 


is  the  Monte  Carlo  estimate  over  the  region  S . ,  and  IL  and  Q,    are  random 
variables  with  distributions  d\i   and  uniform  over  S.,  respectively. 
Neglecting  correlations  between  regions  S.  and  S.,  we  obtain  for  the  total 
variance : 

K 


var^)  ?  S  [A(S  )]"  var(<f>M  ) 


(5.2.20) 


and 


■k.'    h    A  varf£(Sk-8kJ) 

M.    k=l 
J 


var«f>M  )  =  — -         S  var{f(5.,Q.(j))} 


+  2  S   cov{f(§  Q(j)),f(?   Q(J})} 
k<k' 
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Since  some  of  the  variables  in  %,,     are  chosen  by  a  Markov  chain,  the  co- 
variance  term  is  not  zero.   For  simplicity,  however,  we  assume  correlations 
between  terms  in  the  sample  are  negligible.   The  variables  Q,    are  all  of 
the  same  distribution  and,  with  the  exception  of  the  variables  generated 
by  a  Markov  chain,  so  are  ?,  .   The  Markov  chain  variables  approach  the 
distribution  of  dP    after  many  transitions  and  we  will  assume  the 
distributions  to  be  equal.   Hence, 

M. 

E  var{f(^,  ,Q,(j))}  ~  M.  var(f .) 
k=1       ~K  ~K        J      J 


where  var(f.)  is  the  variance  of  f  over  region  S., 

var(f  )  =  ^-y  J  dQ  J  d^f  (?,Q)]2 

j   S. 
J 

"  [A(SJ  J  d2  J  d^  f(|>S)]2  *  (5.2.21) 

Thus,  equation  (5.2.20)  becomes 

K   [A(S  )]2 

var^)  r  E  ^ var(f)  .  (5.2.22) 

j=l     J 

As  an  approximation  to  the  minimization  of  the  actual  variance  of  g^,  then, 

we  choose  M.  so  as  to  minimize  (5.2.22)  subject  to  the  condition  that 

K         J 

E  M.  =  M.   Regarding  the  corresponding  continuum  problem  of  minimizing 
j=l   J 

K 
T  =   E  C./x.  (5.2.23) 


subject  to  the  constraint 


1-1   J   J 


K 

E  x.  =  M  ,  (5.2.24) 
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we  find  that 

/C 
x.  =  M  - 


2  /C. 
i=l   L 


is  a  necessary  condition  for  the  minimization  of  T  subject  to  (5.2.24). 
Evidently,  for  the  integer  case,  if  we  let 


A(C.)  /var(f.) 


M.  =  int  1      M 


,       K 


(5.2.25) 


S  /var(f.)  A(S.) 
i=l       X     X 


then  (5.2.22)  will  be  close  to  its  minimum. 

5.2.3  Multi-Stage  Sampling  Scheme 

We  desire  to  employ  equation  (5.2.25)  for  the  evaluation  of 
n_,,(r,M)  in  equation  (5.2.17).   Unfortunately,  the  variances  (5.2.21) 
are  not  known.   We  therefore  perform  sampling  in  stages.   In  each  pass,  an 
estimate  of  the  variance  in  each  region  S.  is  made,  and  from  this  inform- 
ation a  new  distribution  of  M.  is  assigned  and  the  process  continues,  the 
estimations  of  M.  approaching  (5.2.25).   The  total  variance  will  not  be 
at  a  minimum  for  any  particular  M  but  is  expected  to  be  close  to  a  minimum 
especially  for  large  M. 

This  multi-stage  sampling  procedure  will  now  be  described  in  detail: 
Consider  the  L-stage  sampling  procedure  for  the  estimation  of 


K  1 

Z     A(Sj)A(S~yJ1  dQjd^f?>Q)  (5.2.26) 

J  —  -1-  J    c 


J 
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by 


K 


L      J 


(*) 


(*)  n(XJ) 


g  "   S  A(S  )  —  E   S    f(^  ',Q 
j=l    J  Mj  X=0  k=l     ~k       ~k 


(5.2.27) 


where 


M.  =   £  M 


(X) 


4=0  J 


We  calculate  g  by  the  following  algorithm: 

1)   Choose  the  initial  sample  distribution, 

K<°>,  I     M<°>=M<0>  . 
1  i-1   J 


2)   Begin  sampling  with  the  initial  sample  distribution  and  calculate 


M 


(0) 


j   k=1    ~*   ~* 


and 


0<»  -  ZJ  [f(5k(0),qk<°'J>)]2 


3)   Estimate  the  standard  deviation  for  each  region  by 


a.  =  A(S.) 
J      J 


M. 


1   r(2) 


1/2 


4)   Calculate  the  new  sample  distribution  by: 


M 


(1)  ..  **   'i 

j     L   K 

£  CT 
i-1   : 


5)  Make  the  next  pass  (X  =  1,2,...L): 

M<*> 
j     j   k=1   ^        „k 
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M 


(I) 


g.(2>->g?2>+    i    CfCS^.Q^)]2 


k=l 


k   '^k 


a. 


A(S  ) 


"     l  GU) 

_U'=o  j  ' 


2n 


G. 


Z  M 
i'=0  J 


(X') 


1/2 


M 


<*-l) 


M  A 


L  K 

S  a. 
1=1  ] 


6)   Repeat  5)  until  4  =  L  and  finally  let 


K   G.  A(S.) 


g  -  s 


M, 


a  = 
g 


K 

S 

j-l 


a. 
M. 


1/2 


(5.2.28) 


where  O     is  the  approximate  standard  deviation  of  g. 

O 

The  error  of  the  estimation  of  a     has  many  origins: 

1)  Finiteness  of  M  sample  (variance  of  variance  estimate). 

2)  Lack  of  independence  of  regional  estimates  G.. 

3)  Lack  of  independence  of  path  variables  due  to  the  Markov  chain. 
Hence  a  is  not  to  be  taken  too  seriously  but  serves  as  a  rough  indicator 


of  the  accuracy  of  the  estimate 
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5.2.4  The  "Correlation"  Term,   n  , (r) 

The  procedure  outlined  in  Section  5.2.2  for  the  calculation  of 

1  2 

n  ,  (r)  does  not  yield  satisfactory  results  for   n  ,  (r)  as  defined  in 

equation  (5.2.13).   The  following  procedure  yields  smaller  variance  in  the 

result,  particularly  as  r  becomes  large.   We  may  re-express  (5.2.5b)  in 

the  following  form: 

n  ,  (r)  =  2TTW^(r)  JJ  d(T  dQ   <  (^— -  1)  (e   ZJ+e   iJ)> 


s,    '        Z    W*(r) 


(5.2.29) 


where  we  have  used  the  relations 


\e     \e  ))  =  {{e  )e     >,  i  =  1,2, 

W*(r)  =  <e  Vl2>  . 

We  recall  that  the  region  S'  =  SVS^,  is  in  the  first  (Q  >  0,  Q  >  0) 

U  y       z 

quadrant.   It  will  prove  useful  to  extend  the  integral  to  include  the 
image  with  respect  to  the  Q  -axis.   Let  us  use  the  following  notation: 

S0  S  [%  >   °'Qz  >  °>Qy  +  (Qz  "  I  r)2  ^  V     <"  S0) 
S0  H  KJy  >  °'%   <  °'Qy  +  (Qz  +  2  r)2  ±  V 


S   =  compliment  of  Sn  with  respect  to  the  first  quadrant 

in  the  Q  ,Q  plane 
y   z  r 


S^   =  compliment  of  Sn  with  respect  to  the  fourth  quadrant 


Then,  S1  =   S   and 
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2ns,(r)  =  2TTW  JJ  ...  =  nW   JJ 
since  the  integrand  is  symmetric  with  respect  to  the  Q  -axis.   The  integral 

-3V23 

(5.2.29)  may  be  separated  into  two  equal  integrals  involving  e      and 
e      respectively.   These  two  integrals  may  be  combined  to  form: 

-pv12  _p- 

2n      (r)    =    2TT  W*<r)        JJ  dO*  dQ    <  (^ -l)e        13>    . 

b  Z  —       —      y        Z       W*(r) 

souso 

+ 

The  region  of   integration  excludes   S    .      Over  the  region,    Sn,    the   term, 

-pv13 

e      ,  is  considered  to  be  zero  for  most  paths  so  that  the  integral  will 
be  unaffected  by  the  inclusion  of  Sn.   Hence  we  can  integrate  over  the 
right  half -plane  except  for  the  region  Sn .   Furthermore,  by  choosing  the 
origin  for  Q  to  be  located  at  the  position  of  particle  2,  the  integral 
assumes  the  simpler  form: 

1 


2        /   x        o     tTB/   n    r     ^2,_    r»  d(cos    8) 
n   ,  (r)   =   8tt  W    (r)  J      Q  dQ  J    — ^-z L 


R0  -i 


X  E{(P(Tj(T),r)    -    l)FCn+(T),Q  +   r)|Tl(l)    =   q(l)    =   0} 


where  F  (Tj  (t)  ,Q  +  r)  was  defined  by  equation  (4.2.6),  and  P(T|(T),r)  = 

B     -1 
[W  (r)]   F(T](T),r).   Finally,  applying  the  multiple  integral  approximations, 

we  have 
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!„s,(r)  =  8nU*(r)J  dQ/Mf^Jd, 


Ro  -1 

,  3n     3n    ,  . 

XJ  ^rJ  ^r  (pD   <*•»•*>  -  1} 

X  D+(x,y,u,v,Q  +  R)  +  9  (-5)  (5.2.30) 

n 


where  D  (x,y,u,v,Q  )  was  defined  by  equation  (5.2.10),  and 


Q  =  (0,Q  /l  -  cos  6,  Q  cos  6) 


2 
Equation  (5.2.30)  is  suggestive.   It  implies  that  n„, (r)  is  the 

integral  over  Q  of  the  average  difference  between  D  with  and  without  the 

influence  of  the  probability  weight  P    .   The  variable,  Q,  cannot  be 

considered  as  part  of  the  average,  since  the  integral  over  the  product 

space  which  includes  Q  does  not  exist.   This  is  seen  by  changing  the 

order  of  integration  so  that  Q  is  integrated  before  any  of  the  paths. 

The  resulting  integral  is  infinite  for  an  infinite  number  of  paths,  and 

hence  the  iterated  integral  diverges.   By  Fubini's  theorem,  the  integral 

over  the  product  space  does  not  exist.   However,  the  space  of  cos  9  may  be 

considered  as  part  of  the  product  space,  (cos  9)  XxXyXuXv,  since 

it  is  bounded  and  the  integrand  of  (5.2.30)  is  finite. 

Let  us  now  express  (5.2.30)  as  a  Monte  Carlo  estimate: 


R- 

r 

R. 


2ns,(r)  r  2ns„(r,M)  =  8tt  J  dQGM(Q,r)  (5.2.31a) 


0 
where  ,, 

G  (Q,r)  =  -  I   [D  (x(k),y(k),u(k),v(k),Q(k)  +  r) 
k=1 

-  D+(X(k),y(k),u(k),v(k),Q(k)  -1-  r)]Q2  ,      (5.2.31b) 
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and  the  subscript,  S" ,  indicates  that  we  are  approximating  S'  by  a 
truncation  of  the  upper  limit  on  Q  at  R_.   The  distributions  of  the  random 
variables  x(k),  y(k),  u(k),  v(k)  were  discussed  in  Section  4.3.   The 
random  variable,  Q(k),  is  formed  from 


,Q  /l  -  a 


Q(k)  =  (0,Q  /l  -  a(k)%  Qa(k)) 

where  a(k)  is  uniformly  distribution  on  [-1,1].   xOO  has  distribution, 

du   ,  and  can  be  taken  to  be  identical  to  the  "trial"  path  of  Section  4.3. 
nX 

The  advantage  of  this  procedure  is  that  if  we  take  the  value  of  5  in 
equation  (4.3.20)  equal  to  zero,  the  value  of  xOO  will  equal  the  value 
of  x(k)  for  all  "trial"  paths  which  "succeed."   In  such  cases,  the  summand 
need  not  be  evaluated  since  it  will  be  zero.   It  is  apparent  that  as  r 
becomes  large  (relative  to  X)  the  probability  density,  P    (x,u,r)  -»  1 , 
and  more  and  more  paths  will  succeed.   Hence  relatively  little 
computational  effort  is  expected  for  larger  values  of  r. 

5.2.5  The  "Superposition  Approximation" 

Finally,  we  discuss  an  approximation  which  reduces  the  calculation 
of  rL   (r)  to  an  integral  of  two-particle  density  matrix  elements  only. 
This  approximation  was  used  by  Fosdick  and  Jacobson  [33,34]  and  has  proved 
to  be  useful  for  large  r  and  high  temperatures.   In  this  approximation, 
the  three-particle  term: 

W3  =  <e  ) 

is  replaced  by  the  product  of  the  two-particle  terms: 

_     -pv    -pv    -pv 

•  WW  W+  =  <e   12)<e   23)<e   13>  (5.2.32) 

in  equation  (5.2.1).   The  factorization  is  expected  to  be  valid  providing 
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the  correlations  between  the  terms  are  negligible.   To  see  qualitatively 
how  this  might  come  about  let  us  consider  the  particles  1,  2,  and  3  moving 
in  random  paths  for  the  various  values  of  Q  in  the  integral.   As  the 
temperature  becomes  large,  the  paths  shrink  to  the  points  1,  2,  and  3 
and  the  correlations  become  small  simply  because  the  particles  are  not 
moving  much.   Also  as  r  becomes  large,  particle  3  tends  to  interact  with 
either  particle  1  o_r  particle  2  but  interacts  weakly  with  both  simultan- 
eously.  That  is,  due  to  the  exponential  nature  of  the  Brownian  paths  one 
can  consider  spheres  about  each  point  of  a  certain  radius,  R,  =  0(X),  in 
which  "most"  paths  occur.   Thus  the  chances  that  the  particles  move  out  of 
these  spheres  are  small.   If  the  range  of  the  two-particle  potential  is  r_ 
we  may  consider  concentric  spheres  about  these  spheres  of  radius  rn  +  R, . 

(J     A. 

This  bigger  sphere  is  the  region  within  which  two  particles  can  interact. 
Now  if  r  >  2rn  +  4R,  the  sphere  of  particle  3  cannot  overlap  with  both  that 

of  particle  1  and  of  particle  2  for  most  paths.   Hence  the  important 

-8  (V   +  V   +  V   )  -8(V   -!-  V   ) 

,.  fc.         P^  12     23    13;       .t.      P^  12    23; 
contributions  to  e  are  either  e  or 

-6(V12  +  V13) 
e  and  V- _  is  close  to  zero  for  most  paths  so  that  correlations 

-3v12     -3V.3 

between  e      and  e     ,  i  =  1,2  are  small. 

Replacing  (5.2.32)  in  equation  (5.2.1)  we  have  in  the  "superposition" 
approximation, 

CO         CO 


nP-\r)    r  n   (r)  =  tt  f  dQ2  f  dQ  W(w"'"-1)  (W_-l)  . 

d       sp      J   y  J   z 
Let  us  consider, 


y "      z 

0    -00 


gsP(r)  s  -wfer =  n ?  dQy  ?  dQz  f (Q+)f (Q_)  (5-2-33) 

0      -co 
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where 

+        + 
f(Q')  =   1  -  W"  . 

As  in  Section  5 . 2  we  can  restrict  the  integration  to  the  first  quadrant  of 

2 
the  0   -  Q  plane,  consider  a  "hard  core"  region  SQ  in  which  f  (Q  )  =  1, 

transform  to  the  Q  ,Q  coordinate  system,  and  truncate  the  integral  beyond 

+ 
a  circular  region  of  radius  R,  in  which  f (Q  )  ~  0: 

g   (r)  =  —  JJ  Q+dQ+Q"dQ"  f  (Q+)  +  £-  JJ  Q+dQV  dQ~f  (Q+)f «}") 

SP     r  S0 

(5.2.34) 
=  g.  (r)  +  gq  (r) 

so     si 

where 

SQ  =  {Q+  >  Q":  Q"  <  RQ,Q+  +  Q"  >  r  >  Q+  -  Q_] 

and 

S-l  =  'Q+  >  Q":  Rq  <Q\<  RrQ+  +  Q"  >  r  >  Q+  -  Q-}  . 

Nov,  a  table  of  f (x)  can  be  calculated  by  Monte  Carlo  sampling  as  suggested 

in  Section  4.2.   The  first  integral,  g   (r)  can  be  expressed  as  a  one 

dimensional  quadrature,  and  the  second  can  be  evaluated  by  an  iterated 

quadrature  method.   The  details  are  discussed  in  Appendix  A. 

The  effectiveness  of  the  "superposition"  approximation  is  demonstrated 

in  the  following  table  of  g   (r)  compared  with  path  integral  calculations 

sp 

done  by  Jordan  [9]  not  employing  the  "superposition  approximation,"  for 
T  =  5  K  and  10  K.   Other  comparisons  are  found  in  L33]. 


Table  1 


Comparison  of  n   (r)  Obtained  by  Superposition  Approx- 
imation with  Path  Integral  Calculations  of  rL  (r) 
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T  =   5°K 

T  = 

:   10°K 

r 

n 
sp 

tf}< 

r) 

n 
sp 

r,(1) 
nD 

(r) 

1.0 

.22 

.05 

+ 

.02 

.30 

.10 

+ 

.02 

1.3 

.20 

.01 

+ 

.07 

.04 

-.18 

+ 

.07 

1.5 

-.02 

.14 

+ 

.07 

-.40 

-.42 

+ 

.05 

1.7 

-.12 

.22 

+ 

.06 

-.51 

-.52 

+ 

.04 

1.9 

-.11 

.14 

+ 

.04 

-.45 

-.50 

+ 

.03 

2.1 

-.04 

.12 

+ 

.03 

-.29 

-.27 

+ 

.02 

2.3 

.06 

.06 

+ 

.02 

-.10 

-.10 

+ 

.02 

2.5 

.14 

.16 

+ 

.02 

+  .05 

+  .08 

+ 

.01 

2.8 

.19 

.21 

+ 

.01 

.130 

.132 

+ 

.006 

BO 


5.3  The  Pair  Exchange  Term 

According  to  equations  (2.4.6-7),  the  contribution  to  the  pair 

distribution  function  due  to  the  exchange  of  two  of  the  three  particles  is 

split  into  two  integrals  n    (r)  and  n    (r) .   The  second  involves  the 

1  2 

exchange  of  particles  2  and  3;  and  in  the  Wiener  integral  formulation  a 

(Q~)2 


Q 

factor  of  e       occurs  which  can  serve  as  a  weighting  function  for  the 
integration  over  the  third  particle  as  in  the  cyclic  exchange  term.   The 

np   (r)  integral  however  involves  the  exchange  of  particles  1  and  2  and 

1  _r2/g 

hence  contains  the  factor  e     ,  and  so  the  convergence  of  the  integral 

over  the  third  particle  (Q)  depends  upon  a  cancellation  between  two-  and 

three-particle  terms.   Thus  n    (r)  is  calculated  by  a  method  similar  to 

1 
that  for  il   (r) ,  and  il  \r)  is  calculated  by  a  method  similar  to  that 

for  n    (r) .   It  is  clear  that  il   (r)  should  be  more  difficult  to  cal- 

1 

(1) 
culate  than  il,   (r) .  Moreover  for  the  purpose  of  calculating  the  third 

2 

virial  coefficient  it  is  not  necessary  to  have  both  pair  terms.   For  these 

reasons,  we  will  not  be  concerned  with  the  calculation  of  n^   (r)  . 

1 

We   eliminate  ru,      (r)   by   considering   the  third  virial   coefficient 
1 

term,   a,    =        d  r  (n^      (r)  +  tl,      (r)).      Combining   the   integrals   over  r    , 

12  3 

P 
using   equations  (2.4.6-7)  ,   a,    can  be   expressed 

c£  =  J*  d3r2d3r3{[\9<2   1   3|e        3|l    2   3>   -    \6<2  l|e"      'l    2> 

X    (\b<2   3|e        l\2   3)  +   Xb<3   2|e        2|3   2»] 

+   2[\y<l    3   2|e        3|l    2   3)   -    T<3   2|e        2|2   3> 

X    (Xb<l    2|e        Z\l    2)  +   \b<2   l|e        2|l    2»] 
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ft  "0H9  ft  -^H9 

-  X6<2   l|e        2|l    2)    (Xb<3   l|e        2|3   l) 

ft  -^H9 

-  \°<3    2|e  |2    3)    -    1)}    .  (5.3.1) 

The  integrand  is  of  the  form,  [...]  4-  2[...]  +  ...,  where  the  terms  in 
the  second  bracket  are  the  same  as  those  of  the  first  bracket  except 
r.  <— >r«.*   Since  the  resulting  integrals  over  r  and  r  of  each  bracket 
are  independent  of  r..  ,  r  ,  and  r  ,  the  variables  r..  ,r  ,r  are  dummy 
variables  and  the  integral  of  each  bracket  is  invariant  under  the  trans- 
formation r..  <  >  r  .   Hence  the  two  brackets  have  equal  integrals  and  we 


may  write 


al  =  U  d'r  ^^^  +  ^  d3r2d3r3WE(r12)(WD(r31) 


-  wE(r23)  "  1)  (5.3.2) 


where 


,6,      -PH2 


Wn(r..)  h  \°<r.r.|e   Z|r.r.)  (5.3.3) 

D   lj        l  j  '      '  l  j 


.6,      "PH2 


W_(r..)  =   X (r.r.le     |r.r.)   .  (5.3.4) 

E   ij        j  l1       l  j 


E.g.  the  term 


—21    —13    —23 

RH.  9   "PV1?   "PV^   "K^ 

(2  1  3|e"pH|l  2  3)  =  X"9(e   12  e    23  e   13> 

—23    —31    —21 

-SV    -SV    -SV 

becomes  \~9<e   32  e   21  e   31)  =  <1  3  2  |e"^H|l  2  3)  since  V(r)  =  V(-r) 
and  each  conditional  Brownian  path  has  the  same  distribution. 


B2 


Due  to  the  presence  of  an  exponential  weighting  function  for  the  third 
particle,  the  integral  over  r~  of  the  two-  and  three-particle  density 
matrices  may  be  separated  out  as  follows.   Let 


v 


9  p  ,3    „  _,  -?H3 


(t12)  s  X  J  d  r3  (1  3  2|e   J|l  2  3>  (5.3.5) 


then 


a^  =  3  J  d3r(v(r)  -  a)  +  2a(2b-a)  (5.3.6) 

where 

a  =  J  d3rW£(r)  (5.3.7) 

b  =  J  d3r(l  -  WD(r))  .  (5.3.8) 


In  terms  of  Wiener  integrals  the  three-particle  term  is 

,3   "  R  (Q  } 
V(r)  =  J  dJQ  e  P 

X  E{F[5»Q+;T]]F[Q",-Q";7]"]F[Q+,r;5+]|5(l)  =  q(l)  =  0}. 

(5.3.9) 

Since  the  integral  is  over  all-space,  we  can  integrate  with  respect  to 
Q  =  Q  -  —  r  just  as  well.   Using  equation  (4.2.9)  and  expressing  (5.3.9) 
in  terms  of  spherical  coordinates:  Q  ,  a  -  cos  8,  and  0  we  have  (dropping 
the  superscript  on  Q  ) 

V(r)  -  2n  J  q2dQ  J  da  ."  »  '  Jd,xnJ%nI^!j^! 
0-1  '22 

xexp[-|  S  J  (v\l  +   V^  +   V^)dr]  +0(^)     (5.3.10) 
i=l  o  n 
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where 


V^   h   V(/23    (i-  p(u,T)    +  x      ,    +  t(x   -x        )) 
/n  ~~  ~i    i  111 


+  r  +(  T+f"L    )Q)  (5.3.11a) 


V^  =  V(/2p    (7-  P"(u,v,t)   +  xT        +  t(x'    -   xT      )) 
/j  /n  ~     ~  ~  ~i-i  ~i       ~i-i 


+  Q  -   2(:LVi)Q)  (5.3.11b) 


V13   ~  V(/2P    (7~  £+("'y'T)   +  *i-l  +  r(*i   "   *i-l)} 


+   Q  +  r   -    (I±|li)Q)  (5.3.11c) 


and 


Q  =  (0,Q/l-a2,   Qa)    . 


The  Q   integration  can  be  normalized  by  expressing   the   exponential   factor 

2 
in   terms   of   the  x     distribution  with   three   degrees   of   freedom: 

_   1     _  3  1         2 

T(Q2)    =   2  tt     2  p      2    (Q2)2  e"Q   /P    .  (5.3.12) 


Then 

2  „  J  Q2dQ  J   da  e"^    . . .    .  ^TJ  J  r(Q2)dQ2  f  f»   . . . 
0-1  2  0  -1 

We   next   choose   as   a  weighting   function   for   importance   sampling, 

n      1 
P^n)(x,y,u,v,Q2,a)   =  \~  exp[-  |     2     J   V^2   dT]  (5.3.13) 

P  1=1   0 
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where 


0-1  ll 

n      1 


X  exp[-  |     S     J   V^3   dr]    .  (5.3.14) 

n  i=1  0 


Hence 


and 


J  T(Q2)dQ2    ...    P<n)(x,y,u,Q2,a)   ^  J  dP^n)   =   1 


*<*>  '  372  zp  I  dppn>  «pC"  I  .=,  J  <vu  +  vn'dT] 

2  1=1   0 


+  0(-~)    .  (5.3.15) 


n 


2' 


Now,    it   is   easy  to  verify 


so 


and 


where 


3/2 
Zp  =  ^-3—  a '+  O(^)  (5.3.16) 

\  n 


V(r)   =  a  U(r)   +  0  (^)  (5.3.17) 

n 


a*  =    12  tt  a  J   r2dr[U(r)-l)]  +  2a(2b-a)  (5.3.18) 

0 


n      1 
U(r)    =  J   dP<n)   exp[-  J     E     J    (v[*  +  vJ*)dT]  (5.3.19) 

i=l   0 

M  n      1 

~.  ±     S     exp[-  I     2     J    (Vj~(k)   +  V^(k))dT]  (5.3.20) 

k=l  i=l  0 
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is  calculated  by  Monte  Carlo  sampling  as  described  in  Section  4.3.   In  this 

2 
case  all  random  variables  x(k),  y(k),  u(k),  v(k),  Q  (k) ,  a(k)  are  obtained 

by  importance  sampling  via  a  Markov  process  based  on  dP    . 

5.4   The  Path  Integral  Programs 

Path  integral  computations  were  performed  on  an  IBM  System/360  model 
75  computer  by  means  of  four  FORTRAN  programs  which  we  shall  identify  as 
C,  D  ,  D2,  and  P  corresponding  to  equations  (5.1.9),  (5.2.17),  (5.2.31a) 
and  (5.3.20)  respectively.   These  programs  were  based  on  the  methods 
developed  in  the  previous  sections.   In  all  programs,  the  Markov  chain 
parameter,  6  (see  equation  (4.3.20)),  was  taken  to  be  zero,  and  the  same 
random  numbers  were  used  for  all  values  of  r  for  a  given  run.   Although  this 
tended  to  correlate  the  results  with  respect  to  r,  it  presumably  saved  con- 
siderable computational  time. 

The  initializing  vector  §(0)  for  the  Markov  chain  mentioned  in  Section 
4.3  was  determined  in  the  following  way.   First,  the  arbitrary  value  of  zero 
was  assigned  to  all  path  components;  next  a  subroutine  used  for  determining 
the  next  vector  in  the  chain  was  called  repeatedly  for  a  number  of  times, 
M„,  until  it  could  be  assumed  that:  (a)  a  "trial"  vector  had  "succeeded" 
for  all  values  of  r  (U(§(M  ))  would  then  be  defined  in  the  program);  (b) 
the  probability  distribution  of  §(Mn)  would  be  close  to  P(^(M  ))  (equation 
(4.3.7)).   The  chain  was  then  continued,  being  used  now  in  the  Monte  Carlo 
sampling.   The  initial  process  in  which  the  chain  is  generated  without 
sampling  is  referred  to  as  the  "initial  relaxation"  of  the  chain  and  tends 
to  reduce  the  variance  of  the  result. 

Since  the  states  generated  by  a  realization  of  a  finite  Markov  chain 
tend  to  distribute  about  states  corresponding  to  a  particular  relative 
maximum  of  the  probability  distribution  it  is  possible,  in  the  case  that 
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there  is  more  than  one  relative  maximum,  that  the  state  space  will  not  be 
covered  adequately  by  only  one  Markov  chain.   To  avoid  this  situation,  the 
averages  over  at  least  ten  independent  realizations  of  finite  Markov  chains 
were  taken,  each  chain  being  allowed  an  initial  relaxation  before  sampling 
was  taken.   This  method  has  been  suggested  by  several  authors  [31,35]. 
An  estimate  of  the  standard  deviation  of  the  result  based  on  the 
independent  chains  was  calculated  in  programs  C   and  P.   In  the  case  of 
programs  D.  and  D9 ,  each  independent  realization  corresponded  to  a  stage  of 
the  multi-stage  sampling  procedure  described  in  Section  5.2.3,  and  the 
standard  deviation  was  estimated  according  to  equation  (5.2.28). 
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6.  RESULTS  AND  DISCUSSION 

6.1   Tests  of  the  Path  Integral  Programs 

6.1.1   Harmonic  Oscillator  Tests 

The  programs  C,  D  ,  D  ,  and  P  were  tested  by  temporarily  assuming  the 
two-body  interaction  to  be  a  harmonic  oscillator  rather  than  the  Lennard- 
Jones  potential  (2.1.2).   That  is,  we  let 

V(r)  =  |  r2  .  (6.2.1) 

With  this  assumption,  the  appropriate  three-body  density  matrix  elements 
can  be  expressed  as  products  of  one-dimensional  density  matrix  elements  of 
the  form: 

<X'|exp[-  -£-   (p2  +  m2U)2q2)]|X>  =  p(X,X',3>  (6.2.2) 

z  0 

which  are  solutions  of  the  Bloch  equation, 

*     1   *2    "vX   ? 
(fo  -  ■£-  Z-z  +  -r—  X  )  P(x,x'  ;P)  -  0 

**      2mo  *x 

subject  to  the  initial  condition: 

£im  d(x,x'  53)  -  S(X-X')  (6.2.3) 

3  -  0 

and  the  boundary  conditions 

lira  o(x,x' ;P)  =  o  . 

X  "•  °° 

The  solution  to  the  above  equation  can  be  determined  and  is  given  by 


°(x,x';3)  =  [^  sinh(^)]-1/2 

0 

™0^  2         2 

X  exp{-  —  [coth(u£)(X  +x'    )    -    2XX'    cosech(UUg)  ]} . 

(6.2.4) 
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Therefore  the  pair-correlation  functions  can  be  evaluated  exactly  as  a  check 
against  the  path  integral  calculations.   Some  of  the  results  for  the 
potential  (6.2.1)  are  listed  in  Tables  2  through  8.   Program  D^  was  altered 
slightly  so  as  to  calculate 

I/   d3r3W^(l,2,3)    =  J   d3Q   E{|fV|T|(1)    =   q(l)    =   0} 


rather  than 
2 


n        (t  } 

-f —  =  J   d3Q  e{|  (F+-i)(F'-i)|-n(i)  =  a(i)  =  o} 

W(r)  S' 


to  facilitate  the  analytical  calculation.   The  value  of  the  "hard  core 
radius,"  Rn,  was  set  to  zero.   The  numbers  after  the  +  signs  in  the  tables 
are  the  computed  standard  deviation.   It  will  be  noted  that  in  Table  7 
the  two  independent  calculations  for  n  =  4  agree  to  within  the  estimated 
standard  deviation,  whereas  in  Table  3  for  n  =  8  the  numbers  agree  to  with- 
in two  standard  deviations.   The  reason  could  be  due  to  the  fact  that  the 
standard  deviation  estimates  in  programs  D..  and  D  are  based  on  correlated 
samples  and  are  therefore  inaccurate  whereas  estimates  in  programs  C  and  P 
are  based  on  independent  samples  . 

In  each  table  except  Table  4,  the  effect  of  increasing  the  number, 
(n-1) ,  of  "break  points"  of  the  paths  is  shown  to  be  a  definite  decrease  in 
error  as  expected.   Results  for  largest  n  are  seen  to  agree  with  the  exact 
values  to  within  the  estimated  standard  deviation  with  the  exception  of 
those  from  program  D-.   Errors  in  program  D  have  a  third  source  besides 
the  finiteness  of  n  and  the  variance  of  the  sample:  namely,  the  quadrature 
error  in  the  integration  over  Q.   (A  trapezoidal  rule  was  used  in  this 
integration.)   The  function,  G  (Q,r),  defined  by  equation  (5.2.31b)  is  a 
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Table  2 

Comparison  of  Calulations  of  X~  [1^(1,2)]"   f  d3r  w!?(l,2,3) 
Obtained  from  Program  D.  and  by  Exact  Expression  Assuming 
Harmonic  Oscillator  Potential.  (T  =  0.6,  y   =  .3/T) 


r n  =  4 exact 

0.1  (.136  +  .004)  x  10"1  .132  x  10-1 

0.5  (.129  +  .004)  X  10'1  .129  X  10_1 

1.0  (.114  +  .003)  X  10'1  .118  x  10"1 

2.0  (.85  +  .02)  X  10"2  .86  X  10"2 

3.0  (.47  +  .02)  X  10"2  .50  X  10"2 


Table  3 

Comparison  of  Calculations  of  X  [W  (1,2)]     d  r„  W  (1,2,3) 

Obtained  from  Program  D..  and  by  Exact  Expression  Assuming 

Harmonic  Oscillator  Potential.   (The  two  columns  for  n  =  8 
represent  independent  samples.)   (T  =  0.4,  y  =    .3/T) 

r        n=2  n=4  n=8  n=8         exact 

0.1  (.74+.03)Xl0_3   (.81+.03)X10"3   .802X10"3 

1.0   (.57+.02)xl0_3   (.64t.02)XlO"3   ( .66+.03)xl0"3   ( . 71+.02)xl0_3   .727X10"3 
2.0   (.46+.02)XlO-3   (.42+.02)XlO-3   ( .48+.02) XlO"3   ( .51+.02)xl0"3   .539xl0"3 


V- 


Table  4 

2         B 
Comparison  of  Calculations  of  g   (r)  =  n  , (r)/W  (r) 

Obtained  from  Program  D  with  Exact  Expression  Assuming 

Harmonic  Oscillator  Potential.   (The  two  columns  for 

n  =  2  represent  independent  samples.)   (T  =  0.6,  y   =  ,3/T) 


gD  (r)A3 

r n  =  2 n  =  2 exact 

0       (,38+.03)Xl0"2     (.24±.03)xl(f2      .436xl0'2 
2.5      (.11-K1)X10~2      (.21+.07)Xl0"2      .285X10"2 
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Table  6 

Comparison  of  Results  of  Program  C  with  Exact  Values 
Assuming  Harmonic  Oscillator  Potential  (T  =  0.2,  y   =  .1/T) 
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n  =  4 


2-1/2gc(r)/X3 

exact 


Zc(r) 


n  =  4 


exact 


0. 

(.17+.02)X10 

1. 

(.13+.01)X10 

2. 

(.73+.06)Xl0 

-4 


-5 


.180X10 


-4 


.146X10 


.789X10 


-5 


(.162+.004)X10 
(.149+.004)X10 


-1 


(.116+.003)X10 


.0165 
.0152 

.0120 
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Table  8 


-1 


Comparison  of  Calculations  of  U(r)  and  a  v(r)  Obtained 
from  Monte  Carlo  Sample  and  Exact  Expression  Respectively 

Assuming  Harmonic  Oscillator  Potential.   (The  two  columns 
for  n  =  4  represent  independent  samples.)   (T  =  0.6 

and  Y  =  0.3/T) 


U(r) 


Monte  Carlo  Calculation  (P) 
n  =  2 


exact 


0.0 
0.5 
1.0 
1.5 
2.0 


(.52+.02)xlO 

(.46+.02)XlO' 

(.32+.02)XlO" 

(.17+.01)X10 

(.73+.05)Xl0" 


-1 


-1 


53x10 


-1 


.47X10 


-1 


33X10 


-1 


.18x10 


-] 


.81x10 


-2 


150+.002 


1533 
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fairly  rapidly  varying  function  with  almost  equal  positive  and  negative 
contributions.   If  the  peaks  are  sharper  than  the  quadrature  spacing, 
considerable  error  could  arise.   Since  the  integral  is  the  result  of 
almost  equal  positive  and  negative  contributions,  this  error  as  well  as 
the  standard  deviation  can  be  quite  significant.   The  large  error  in 
column  2  of  Table  4  (r  =  0)  is  due  to  the  fact  that  G  (Q,r)  was  not 
evaluated  at  a  point  near  enough  to  its  positive  peak. 

6.1.2  Moderate  Temperature  Tests 

The  availability  of  independent  calculations  of  n    (r)  at  moderate 
temperatures  [9]  permitted  a  further  check  on  programs  D,  and  D  using  the 
Lennard-Jones  potential  (2.1.2).   Table  9  shows  a  comparison  at  T  =  20  K 
with  calculations  of  Jordan  which  also  used  the  path  integral  method. 
Let 

gD(r)  h  n^)(r)/W^(r)  =  g])    (r)  +  gD  (r)  (6.2.5a) 

where 

g_  (r)  =  (n-  (r)  -I-  Ln_ ,  (r) ) /W* (r)  (6.2.5b) 

Dl      bo  b  l 

gD  (r)  =  2ns, (r)/W*(r)  .  (6.2.5c) 

It  is  seen  that  g   (r)  is  significant  at  r  =  1  and  is  necessary  for  agree- 

2 
ment  with  Jordan's  calculations.   For  the  other  values  of  r,  g   (r)  is  of 

the  order  of  its  variance  and  is  too  insignificant  and  uncertain  to  be 

included.''   Included  for  comparison  are  the  values  obtained  via  the 

superposition  approximation  described  in  Section  5.2.5. 


The  function,  g   (r)  differs  from  Jordan's  g  (r)  in  that  g   (r)  does  not 
2       ±  6  2 

include  the  regions  S_  in  its  integration. 


H 


Table  9 

Results  of  Programs  D..  and  D  Assuming  Lennard-Jones 

Potential  at  T  =  20°K  Compared  with  Calculations  of 
Jordan  and  Values  Obtained  from  "Superposition 
Approximation."  (n  =  2,  R..  =  0.5) 

\  r2 

r     gD  (r)/X3    gD  (r)/\3     gD(r)/\3   gD(r)/X3  (Jordan)   gsp(r)/\3 

1.0   1.048+.066   -.196+. 054    .85+. 08     .78+. 08  1.56 

1.2    .446+. 063    .031+. 027    .45+. 07     .30+. 08  .48 

1.4   -.329+. 057    .020+. 021   -.33+. 06    -.31+. 07  -.275 
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6.2      Low  Temperature  Results 

A  complete  set  of  tables  of  numerical  results  for  the  temperatures 

1.0  K  and  1.7  K  are  to  be  found  in  Appendix  C.   The  functions,  g   (r) , 

1 
g   (r) ,  U(r),  and  gr(r)  were  obtained  directly  from  the  path  integral 

D2  L 

programs  D.  ,  D  ,  P,  and  C  respectively.   Values  of  the  two  particle  Slater 

B         E 
sum's  direct  and  exchange  parts,  W  (r)  and  W  (r)  necessary  for  the  cal- 
culation of  the  correlation  functions,  n^   (r)  and  n^   (r),  were  obtained 
by  additional  path  integral  programs  for  T  =  1.7;  and  in  the  case  of 
T  =  1.0,  accurate  phase  shift  calculations  due  to  Larsen,  Witte,  and 
Kilpatrick  [36]  were  used.   Integrals  of  the  correlation  correction  terms 
for  the  third  virial  coefficient  were  computed  by  means  of  the  trapezoidal 
quadrature  rule.   The  constants,  a  and  b,  necessary  in  the  calculation  of 
0_  (5.3.18)  were  obtained  from  accurate  phase  shift  calculations  of  the 
direct  and  exchange  contributions  of  the  second  virial  coefficient  due  to 
Boyd  and  Larsen  [12]. 

Values  of  g   (r)  (6.2.5c)  are  seen  in  Tables  10  and  11  to  be  distributed 
somewhat  equally  about  the  value  zero  for  both  T  =  1  and  1.7,  and  most  values 
seem  to  be  of  the  order  of  the  estimated  standard  deviations .   Their  con- 
tribution to  the  total,  g  (r),  therefore,  will  be  neglected  so  that  the 
adopted  values  of  gn(r)  are  taken  to  be  identical  to  those  of  g   (r) ,  the 
only  effect  of  g   (r)  being  to  raise  the  estimated  standard  deviation  of 
gD(r). 

Comparison  of  the  adopted  values  of  gn(r)  with  those  obtained  via 

the  superposition  approximation  g   (r)  in  Figure  7  shows  agreement  to 

sp 

within  the  estimated  standard  deviation  for  r  >  1.5  at  both  temperatures. 
There  seems  to  be  a  consistent  difference  between  the  two  functions  for 
T  =  1.0.   This  is  perhaps  not  entirely  due  to  the  correlation  effect 


J.1 
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Figure  7.   Three-Particle  Boltzmann  Contribution  to  the 
Pair  Correlation  Function,  X  8n^r^» 
Compared  with  X  g   (r) 
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mentioned  in  5.4  since  the  calculations  were  in  fact  performed  in  piecemeal 
fashion.   The  complete  function,  g  (r) ,  used  in  the  computation  of  CL  was 
obtained  by  joining   in  a  composite  table,  linearly  interpolated  values 
of  gn(r)  obtained  from  the  path  integral  programs  with  the  values  of 
g   (r)  obtained  by  superposition  approximation  at  the  points:  r  =  4.5  for 
T  =  1.0  and  r  =  2.3  for  T  =  1.7. 

Comparisons  among  the  various  contributions  to  the  total  three- 
particle  correction  term,  n    (r)  ~  n^   (r)  +  n    (r)  +  n    (r)  in 
Figures  8  and  9  show  that  the  direct  term  dominates  for  all  r.   Further- 
more, it  is  seen  that  n   (r)  is  smaller  than  n  (r)  for  intermediate  r  but 

F2  L 

has  a  longer  tail  which  accounts  for  the  fact  that  the  integral,  ol  ,  is 

much  larger  than  a      (see  Table  16)  .   Essentially  two  factors  account  for 
the  smallness  of  the  exchange  terms: 

(a)  the  kinetic  energy  required  for  the  motion  of  particles  in  an 
exchange  (giving  rise  to  an  exponential  factor)  . 

(b)  the  repulsive  core  of  the  interactions,  which  limits  the 
possible  paths  in  which  particles  may  undergo  exchanges. 

Cause  (b)  should  be  especially  important  for  pair  exchanges  (particularly 

as  the  temperature  becomes  higher  and  the  paths  approach  a  straight  line) 

and  accounts  for  the  relative  smallness  of  n   (r)  at  intermediate  r. 

2 
The  exponential  term  of  n  (r)  for  large  r  is  due  primarily  to  (a)  and 

accounts  for  the  shorter  tail  of  n  (r) .   The  term  n   (r) ,  which  has  not 

1 
been  calculated,  is  influenced  strongly  by  both  (a)  and  (b)  and  hence 

should  be  smaller  than  either  n  (r)  or  n   (r)  for  most  r.   Not  inconsistent 

C        P2 
with  this  claim  is  the  fact  that  the  ratios  of  the  integrals  CL,   :  CL,   :  (X 

are  7:138:53  for  T  =  1.7  and  130:1600:1080  for  T  =  1.0. 
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Figure  8.   Direct  and  Exchange  Contributions  to 


n^1)(r)/X3  at  1°K 
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Figure  9.   Direct  and  Exchange  Contributions  to 
n^1}(r)/\3  at  1 .  7°K 
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The  three-particle  term,  n   (r),  is  plotted  in  Figure  10  for 
T  =  1°K  and  T  =  1.7°K  and  compared  with  values  at  T  =  5°K  and  T  =  10°K 

[9].   It  is  observed  that  the  function  becomes  rapidly  large  with  de- 

o        o  * 
creasing  temperature,  and  the  relative  minimum  exhibited  at  5  K  and  10  K 

at  r  ~  1.8  becomes  a  slight  inflection  near  the  relative  maximum  at  1.7  K 
and  1  K.  The  relative  minimum  has  been  interpreted  as  the  valley  between 
the  two  peaks  which  correspond  to  preferred  positions  of  three  particles 
arranged  linearly  to  be  within  each  other's  potential  well  [9],   (Particles 
located  at  r  =  0,1.2,2.4)   As  the  temperature  becomes  lowered,  we  expect 
the  distribution  of  the  positions  of  the  particles  to  become  more  spread 
out  (X  becomes  larger)  so  that  the  region  between  the  preferred  positions 
becomes  an  overlap  region  perhaps  accounting  for  the  appearance  of  the  peak 
at  lower  temperatures.   One  would  expect  the  peak  in  n    (r)  to  correspond, 
on  the  other  hand,  to  a  triangular  arrangement  of  three  particles  since 
it  would  require  the  least  kinetic  energy  for  the  cyclic  exchange  of  the 
particles.   This  arrangement  could  also  account  for  the  peak  in  ru   (r) 
at  the  low  temperatures . 

The  computed  values  of  the  third  virial  coefficient,  C(T),  obtained 
via  equation  (2.3.5)  are  in  qualitative  agreement  with  experimental  measure- 
ments as  shown  in  Figure  11,  i.e.  they  are  negative  and  decrease  with  de- 
creasing temperature.   The  experimental  data  however  are  rather  uncertain 
and  marginal  at  best.   The  computed  value  at  1.7  K  agrees  with  a  smooth 
extrapolation  from  Keller's  data  [37,38],  although  values  extrapolated 
from  Keesom's  data  [39]  seem  to  disagree  with  the  calculation  by  more  than 

an  order  of  magnitude.   The  large  relative  error  in  this  calculation  is 

2 
due  to  a  large  cancellation  from  0L  -  ccn. 


It  has  been  shown  [9]  that  n2   (r)  exhibits  this  dip  at  r  ~  1.8  for  the 
temperature  range  5°K  <  T  <  273°K. 
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(1)     3 
Figure  10.   n    (r)/X  as  a  Function  of  Temperature 
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Figure  11.   Computed  and  Measured  Values  of  the  Third 
Virial  Coefficient,  C(T),  as  a  Function  of 
Temperature 
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A  theoretical  calculation  due  to  Larsen  [7]  based  on  the  binary 
collision  expansion  method  of  Lee  and  Yang  [5]  yielded  the  result, 
C  (1.7)  =  -2.2  X  10   cm  /mole  which  deviates  by  over  an  order  of  magnitude 
from  the  result  of  the  present  path  integral  calculation,  C     (1.7)  = 
-8.8  x  10  cm  /mole  .   In  addition  to  using  a  low  temperature  expansion, 
Larsen' s  calculations  were  based  on  the  following  assumptions: 

(1)  that  the  two-body  interaction  potential  is  given  by  a  square 
well  with  a  finite  repulsive  core, 

(2)  that  there  exists  a  three-body  bound  state,  and 

(3)  that  the  binding  energy  could  be  approximated  by  neglecting  the 
repulsive  core. 

Assumption  (1)  can  account  only  for  a  minor  difference  between  C  (1.7)  and 

1j 

C     (1.7).   Assumption  (2)  is  quite  important,  requiring  the  addition  of  a 

tern  proportional  to  e     where  E  is  the  binding  energy;  and  it  is  possible 

B 

that  the  error  resulting  from  assumption  (3)  is  considerable. 

6.3   Conclusions 

In  summary,  the  results  of  the  previous  section  suggest  that  the 
exchange  contributions  to  the  pair  distribution  function  for  He 
gas  are  almost  negligible  at  as  low  a  temperature  as  1.7  K  but  are 
significant  at  1  K .   The  third  virial  coefficient,  however,  due  to  a  large 
cancellation,  is  sensitive  to  three  particle  exchange  effects  at  T  =  1.7  K 
as  well  as  at  1  K.   It  has  also  been  shown  that  three-particle  contributions 
to  the  pair  distribution  function  not  involving  exchanges  can  be  reasonably 
approximated  in  terms  of  two-particle  diagonal  conf igurational  density 
matrix  elements  through  the  "superposition  approximation"  for  temperatures 
as  low  as  1  K. 
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The  success  of  the  superposition  approximation  suggests  the 
feasibility  of  a  so-called  "expanded  superposition  approximation"  of  the 
form: 

<r»  |e     |r  >  =  J  d  xx  ...  J  d  x^  n   <x.+1|e     |x.)  ; 

i=0 

N    N   N     ,N 
*0  =  I    '  *n  =  I 

~T  ~a     n-1  N 

=  r  d3NXl  ...  f  d3Nx  .  n  n 

J     l  J     n"1  1=0  j«c 

x^iik)i«    2     i*--J'k)>  +  E„ 

where 


and 

x.(j'k)  =  (X.   .,  X.  .) 

~i        ~i,j'  ~i,k 

N 
is  a  6- tuple  consisting  of  the  components  of  x.  corresponding  to 

particles  j  and  k.   By  expressing  the  density  matrix  elements  in  terms  of 

Wiener  integrals  and  expanding  the  exponential  functionals  in  powers  of 

the  paths,  one  can  show  that  the  error  term,  E  ,  of  this  approximation 

is  0(l/n2). 

There  seems  to  be  no  obstacle  to  the  attainment  of  further  path 

integral  calculations  of  the  third  virial  coefficient  at  lower  temperatures 

except  an  increase  in  the  necessary  computing  time.   It  might  be  possible 

that,  even  below  1  K,  a  certain  reduction  of  effort  could  be  attained 

through  approximation  of  the  direct  term  by  the  "superposition 
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approximation."   It  should  be  relatively  straightforward  to  extend  the 
path  integral  approach  to  the  case  of  non-additive  potentials;  and  in 
principle,  the  path  integral  formalism  could  be  employed  for  calculations 
involving  four  or  more  particles,  although  computing  time  would  be  a  severe 
limiting  factor.   With  the  use  of  coming  parallel  processing  computers, 
this  limitation  can  probably  be  considerably  eased. 


APPENDIX  A 


QUADRATURE  FOR  g_  (r) ,  gQ  (r) 

so    bi 


We  have,  by  equation  (5.2.31), 
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gSP 


(r)  r  gq  (r)  +  gq  (r) 


where 


gs  (r)  =  ^U   Q+dQ+Q"dQ'f(Q+) 


0 


gs  (r)  =  &  JJ  Q+dQ+Q"dQ"f(Q+)f(Q*) 
1        S, 


and 


S„  = 


.+ 


.+ 


[Q  >  Q"  :  Q"  <  RQ,  Q  -!-  Q"  >  r  >  Q  -  Q-} 


=  {Q+  >  Q"  :  RQ  <  Q"  <  RL,  Q+  +  Q_  >  r  >  Q+  -  Q-} 


We  consider  first  g   (r) .   It  can  be  verified  that 

r+R0 
r   |2  J"   xdx(RQ-(x-r)2)f(x)  +  2tt[|  Ro-Rq(|)  +  |(f)3] 


gs  (D  =  < 


R 


r+R, 


2TT 

r 


xdx(RQ-(x-r)2)f(x) 


if  2  <  R0 


if  1  *  Ro 


(A.l) 


r-R, 


where  use  is  made  of  the  fact  that  f (x)  =  1  for  x  <  R_.   The  integrals  in 

equation  (A.l)  can  be  evaluated  by  a  modified  trapezoidal  rule  making 

2      2 
use  of  the  weighting  function  x(R_-(x-r)  ).   Consider  then,  the  integral 

over  the  interval  [x, ,x  ]: 


I(XVX2)   -  J   xdx(R2-(x-r)2)f(x) 


and  let  us  Linearly  interpolate  f(x)  within  this  interval,  i.e., 
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x   -  x  \  f  x  -  x, 

f  (x)  r   |  — f  (x.)  + 

1  x„  -  x„        1        X„  -  X, 


'2    1 


'2    1 


f(x2) 


Integrating,  we  obtain, 


I(x1,x2)  = 


R0  "  "' 

^(x^xp  +|02(Xl,x2) 


20  03(X1'X2} 


X  (f(Xl)  -  f(x2)) 


+ 


R2  -  r2 

10 1_  ,  2   2,    2   ,  3   3,    1,4  4. 

- (x2~x1)  +  -  r(x2-x1)  -  4^x2~xl') 


f(x2) 


where 


2     2 
J1(x1,x  )  =  x1x2  +  x2  -  2xL 


2        2    3     3 
02(xisx  )  =  x^  +  xxx2  +  x2  -  3xL 


*    r  ^-3   ^22j     3^4 

3   l'X2^  =  X1X2   X1X2   X1X2   X2 


-  4x, 


Then  assuming  a  table  x.,  f(x.)  for  i  =  1 , .  .  . ,N  is  known,  we  can  express 

the  integral  over  an  arbitrary  interval  [a,b]  as: 

b 
I(a,b)  =  J  xdx(R2  -  (x-r)2)f(x) 


V1 

E   I(x.  ,x   )  +  I(a,x.  )  +  I(x.  ,b) 
i=j  Jl       J2 
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where  x.   ,  <  a  <  x.   and  x.   <  b  <  x.  ,,  and  f(a)  and  f(b)  can  be  deter- 

mined  by  linear  interpolation. 

Next  we  consider  g   (r) .   There  are  three  cases. 

bl 

Case  I:  ^  <  RQ 

R 

1  0  +R 

Sc;  (r)  =  ^  J  Q"dQ"f(Q")  J   Q+dQ+f(Q+)  - 

1  R  Q 

0  * 


Case  II:  R0  <  f  <  RL 


.        r/2  r-hQ" 

<r>  -  f1  J   Q'dQ'f (Q")  J"   QdQ+f(Q+) 

R0  r-Q" 

4tt  n1  -  -   -  Vr  +  +   + 
+  f1    J  Q  dQ  f(Q  )  J   Q+dQ+f(Q  ) 

r/2  q- 


g 


s. 


Case  III:   R.  <  | 


Rl  r+Q~ 

gs  (r)-~J  Q'dQ'f (Q")  J   Q+dQ+f(Q+) 

1         R0  r-Q" 


The  above  integrals  are  of  the  form 
r 


2         x+r 
^J"  xdxf(x)  J*   ydyf(y)  .  (A.  2) 

r,         a  (x) 


Define  accordingly, 


i|f(x1}x2)  =  [x  f(x2)  +  Xj^f  (x1)](x2  -  x^ 


x2 


which  is  the  trapezoidal  rule  approximation  to  2  J   xf  (x)dx.   Then  (A. 2) 

Xl 
can  be  approximated  by: 


Ill 


J2"1 


J  {      S       [(g(x.)Y(x.)   +  8(x1+1)Y(xi+1))(xi+1-xt)] 


i=J 


+    (g(r    )Y(r   )   +  g(x      )Y(x      ))(r,-x      ) 

^  ^  Jo  J  9  ^         Jo 


+    (g(x.    )Y(x.    )   +  g(r    )Y(r    ))(x.    -r    )} 
Jl  Jl  Jl 


where 


g(x)    =  xf  (x)    , 


XJ1"1-rl<XJ1    '      XJ2"r2<XJ2+l    ' 

k2-l 

Y(x)    s  iKx.,x        )   +   \|r(a(x),x^    )   +   tCx,     ,x  +  r) 


and 


"k^l   ^a(x)    <Xkx    '      \2^X+   r   <Xk2+l    ' 
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APPENDIX  B 
STRATIFIED  SAMPLING  -  ANALYSIS  OF  REGIONS 


We  consider  the  sets  S.,  i  =  1,...,7  defined  in  Section  5.2.2  by 

s.  =  st  n  s"  s.  =  st  n  s" 

111  431 

52  "  4  "  Si  S5  "  S3  "  S2 

53  "  S2  "  S2  S6  "  S3  °  S3 

S?  =  7   0  (S"  U  S"  U  s-) 


where 


and 


t   ,  +  -  + 

Si  =  [Q  ,Q   :  Rt_x  <  Q  <  Q  <  R.}  n  S 


S  =  {Q+,Q~  :  Q+  +  Q"  >  r  >  Q+  -  Q";  Q+  >  Q"}  . 


Now  the  sets,  S. ,  can  be  expressed  in  terms  of  two  parameterized  sets  as 


follows 


Sl  ~  SI(R0'R1)  S4  "  SII(R0»R1;R2'R3) 

$2  =  Sjj  (Rq>R^'jR-^)R2^    ^5  =  ^ii  (R^>R2'R2'R3^ 
S3=  SjC^.Rj)  S6=  Sl(R2,R3) 

Sy  =  S   (R  ,R  ;R  ,r  +  R«) 


where 


if  Rx  >  R2  or  R2  <  | 


SI(R1,R2)  = 


[Q0  <  Q"  <  Q3;  Q2  <  Q+  <  Q31  n  S  otherwise 


where 


and 


%  =  h>   Q2  = 


Q3  =  R2 


113 


!•*!<! 


R1  ,  otherwise 


0  if  RQ  >  Rx  or  R  >  R  or  R^R,  >  r  or  R3+R,  <  r 


SII(R0'R1;R2'R3)  " 


+ 


[Q0  <  Q"  <  Q1;  Q2  <  Q  <  Q3]  n  s  o 


therwise 


where 


r  -  R3,  RQ  +  R3  <  r 


%=\    R2  "  r'  R2  "R0>r 


R_,  otherwise 


Qi-  \ 


Q2  = 


*3  = 


r  -  R^  Rx  h  R2  <  r 


R  ,  otherwise 


r  +  R1S  R_  -  R.  >  r 


R„,  otherwise 


Th 


+.2 


e  normalization  factors,  a{S.}  =  J   d  (Q  )   d(Q),  can  be  calculated  by 


determination  of  A{s  (R  ,R  )}  and  a{s   (R^LjR  R  ) 
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As  an  example  of  how  the  A's  can  be  determined  we  consider  S  (R.  ,R  ) 
as  shown  in  Figure  12.   The  shaded  area,  S  (R- ,R  ),  can  be  regarded  as  the 
intersection  of  the  rectangular  region  (S)  and  the  isoceles  right  triangle 
in  Q  ,Q  space.   All  possible  "areas"  can  be  built  up  by  a  superposition 
of  "elementary  areas,"  A{cr.},  e.g.  in  the  figure, 
A{SI(R1,R2)}  =  Ata^r-R^}  +  A{a2(R1,R2)} 


where 


Figure  12.   The  Regions,  S,  S  and  a. 


(S  s  (Q+)2,  Tl  ■  (Q")2) 


,2 


aCqjCR')}  =  J   d|     J     dT|  =  rR'Z(|  R'-r)  -  \ 


1   4 
r 


,r.2     ,  .1/2.2 

2 
R2       F 


?       R   -  (r-R  ) 

A{a2(R1,R2)}  ■    J   2  d§  |  dTl ^ — « 1 R2(R2   (r.Ri)2 

(r-Rx)     R2 


1 v  2 
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The    calculation   of   all    possible    regions    S    (R     R   )    and    S      (R     R    ;R     R   ) 
is   achieved   by    the   determination  of   a    total    of  nine   "elementary   areas." 
They   are    (including   the    two   just   defined), 

f  K/2>2        8        3 

Ata,(R)}  -  J     dT|  J  =  f  rR 

0  Cr-lf/2>2 

2  2 

(R2-rT  R2 

A{a    (R     R    )1   =  J"  dTl  J*  d|   =    (R2-r2)[(R   -r)2-R2] 

2  1/2   2 

R^  (r+Tl17    ) 


|[(R2-r)4   -  Bj]   -  |  r[(R2-r)3   -   rJ] 


2  2 

R2  R2 

A{a5(R15R2)}   -J     dTl  J     d?   =  \    (R2   -   R^)2 

R2    1 


r//2    R/2-y 
A{CT  (R)}  =  2  j    dy  J     dx(x  -yZ) 


://2 


+   1  -   1 

where  Q   =  —  (x  +  y)    Q   =  —  (x  -  y) 

/2  /2 


=  Afj^R)} 


2         2 
R2       R3 


A{a7(R15R2,R3)}   =  J     dTl  J 


2  1/2   2 

R^  (r-Tf'V 


=    (R2    -    r2)(R2    -   R2)   +  |  r(R^    -    R^)    -   |    (R*   -   rJ) 
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2     2 
R2    R4 
A{a8(R1,R2;R3,R4)}  -J  dT\  J  d§  =  (R*  -  R^CR*  -  R*) 

Rl    R3 

2  2 

R3        Rl 
A{a9(R19R2JR3)}  =  J*  d§      J"    =  A{a7(R2,R3,R1)}  . 

R*    (?1/2-r)2 

Flow  charts  for  the  calculation  of  A{s  (R  ,R  )}  and 
A{S   (R0,R, ;R_,R~)}  in  terms  of  A{a.}  are  exhibited  in  Figures  13  and  14 
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APPENDIX   C 

TABLES   OF   NUMERICAL   RESULTS   FOR  TEMPERATURES 
1.0°K   AND   1.7°K 


Table    10 

Pair   Distribution  Function:   Three  Particle   Direct  Con- 
tribution for  T  =   1.0°K 


n=19,    RQ=0.75 
r  gD    (r)/\3  gD    (r)/X3        n^l)(v)/\3        g      <r)/\3     ng    (r)/\3 


0.9 

.690+. 06 

-.05+. 06 

.114+. 01 

1.0 

.594+. 04 

-.08+. 07 

.354+. 03 

.8596 

.5118 

1.5 

.570+. 06 

+.04+. 06 

1.28+.1 

.6104 

1.3753 

2.0 

.510+. 06 

.95+.1 

.4757 

.8907 

2.5 

.464+. 06 

.664+. 08 

.4254 

.6085 

3.0 

.398+. 03 

.471+. 04 

.3657 

.4322 

3.5 

.305+. 02 

.326+. 02 

.2726 

.2915 

4.0 

.189+. 01 

.194+. 01 

.1774 

.1819 

4.5 

.116+. 01 

.117+.  01 

.1036 

.1046 

Table  11 

Pair  Distribution  Function:  Three  Particle  Direct  Con- 
tribution for  T  =  1.7°K 
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n=10, 

R0=0.7 

r 

gD    (r)AS 
1 

gD    (r)A3 

n^CD/X3 

ssp(r)/x3 

n      (r)/X3 
sp 

0.9 

.411+. 04 

.03+. 07 

.066+. 01 

1.0 

.421+. 04 

-.1+.07 

.2 25+. 04 

.755 

.404 

1.45 

.336+. 04 

.03+. 04 

.651+. 08 

.466 

.902 

2.0 

.304+. 03 

.07+. 04 

.479+. 06 

.320 

.504 

2.6 

.322+. 02 

-.07+. 04 

.388+. 03 

.321 

.387 

3.0 

.320+. 02 

.02+. 02 

.346+. 02 

.294 

.317 

3.4 

.243+. 009 

-.03+. 01 

.251+. 01 

.230 

.238 

3.6 

.186+. 01 

.190+. 01 

.193 

.197 

3.8 

.162+. 01 

.164+. 01 

.157 

.159 

4.0 

.1394+. 008 

.1406+. 008 

.125 

.126 

5.0 

.0259+. 002 

.0260+. 002 

.030 

.030 
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Table  12 

Pair  Distribution  Function:  Three  Particle  Cyclic 
Exchange  Contribution  for  T  =  1.0°K 


T=l 

.0°K,  n=19 

r 

gc(D 

Zc(r) 

nc   (r)/X 

0.8 

.072+. 02 

.09+. 02 

.0035+. 002 

0.9 

.431+. 07 

.362+. 05 

.079+. 02 

1.0 

.83+. 14 

.645+. 03 

.251+. 06 

1.1 

.97+.1 

.898+. 02 

.378+. 05 

1.2 

1 

.36+.1 

1.059+.04 

.568+. 07 

1.4 

1 

.54+. 16 

1. 288+. 05 

.634+. 09 

1.6 

1 

.81+. 1 

1.271+.02 

.577+. 04 

1.8 

1 

.83+. 09 

1.  262+. 02 

.440+. 03 

2.0 

1 

.926+. 06 

1.231+.02 

.331+. 02 

2.2 

1 

.93+. 07 

1.212+.02 

.232+. 01 

2.4 

1 

.925+. 07 

1.204+.02 

.159+. 009 

2.6 

1 

.926+. 06 

1.200+.02 

.106+. 005 

3.0 

1 

.91+. 07 

1.186+.03 

.042+. 003 

3.5 

1 

.92+. 1 

1.159+.02 

.0110+. 001 

6.0 

1 

.355+. 04 

0. 996+. 06 

(.44+.04)xl0-6 

8.0 

1 

.01+. 02 

0. 928+. 01 

(.36+.01)XlO-11 
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Table  13 

Pair  Distribution  Function:  Three-Particle  Cyclic  Exchange 
Contribution  for  T  =  1 . 7°K 


T=1.7°K,  n=10 

r 

gc(r) 

Zc(r) 

nc(r)/X3 

0.85 

.100+. 01 7 

.1049+. 005 

.0045+. 001 

0.90 

.158+. 03 

.1773+. 007 

.0114+. 002 

0.95 

.233+. 03 

.248+. 009 

.0220+. 004 

1.00 

.380+. 05 

.339+. 006 

.0457+. 007 

1.05 

.391+. 05 

.393+. 02 

.0509+. 009 

1.10 

.472+. 06 

.444+. 02 

.064+. 01 

1.15 

.535+. 05 

.484+. 02 

.074+. 01 

1.20 

.740+. 05 

.537+. 02 

.104+. 01 

1.30 

.741+. 04 

.572+. 015 

.0935+. 007 

1.40 

.924+. 06 

.651+. 015 

.110+. 009 

1.60 

.934+. 05 

.694+. 015 

.0785+. 006 

1.80 

1.026+.04 

.714+. 02 

.0556+. 004 

2.00 

1.125+.07 

.722+. 01 

.0364+. 003 

2.2 

1.203+.04 

.731+. 01 

.0221+. 001 

2.4 

1. 25  2+.  08 

.725+. 01 

.0121+. 001 

2.6 

1.307+.03 

.731+. 01 

(.640+.02)xl0~2 

2.8 

1.317+.05 

.742+. 01 

(.311+.02)X10-2 

3.0 

1.351+.05 

.743+. 01 

(.144+.008)X10"2 

3.2 

1. 345+. 03 

.755+. 01 

(.618+.02)XlO~3 

3.4 

1.340+.03 

.753+. 01 

(.247+.01)XlO~3 

3.6 

1. 336+. 03 

.751+. 01 

(.938+.05)xl0~4 

3.8 

1. 274+. 03 

.747+. 01 

(.321+.01)X10~4 

8.0 

0. 768+. 01 

.692+. 01 

(.263+.009)Xl0'19 

Table  14 

Pair  Distribution  Function:  Three-Particle  Pair  Exchange 
Contribution  for  T  =  1.0°K 
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T=l 

.0°K,  n= 

=  19 

r 

U(r) 

r^  (r)/\3 
2 

r 

U(r) 

n   (r)/X3 
2 

0.8 

.009+. 003 

-.002+. 001 

4.2 

1.0691+.009 

.0270+. 005 

0.9 

.273+. 007 

.039+. 039 

4.4 

1.0407+.005 

.0150+. 002 

1.0 

.89+. 13 

.095+. 07 

4.8 

1.0195+.002 

(.71+.1)X10"2 

1.1 

1.54+.27 

.08+. 14 

5.2 

1.0114+.001 

(.438+.07)xl0_2 

1.2 

2. 49+. 14 

.279+. 07 

5.6 

1.00543+.0006 

(.188+.03)X10-2 

1.4 

2. 91+. 12 

.205+. 06 

6.0 

1.00298+.0004 

(.95+.2)XlO-3 

1.6 

2.81+.1 

.163+. 05 

6.4 

1.00188 

.587X10"3 

1.8 

2.54+.1 

.137+. 06 

6.8 

1.00117 

.345X10"3 

2.2 

2. 18+.  2 

.22+.  12 

7.2 

1.00076 

.211X10"3 

2.6 

1.75+.06 

.186+. 03 

7.6 

1.00049 

.125X10"3 

3.0 

1.52+.05 

.175+. 03 

8.0 

1.00031 

.687X10"4 

3.4 

1. 295+. 02 

.110+. 01 

8.5 

1.00018 

.271X10'4 

3.8 

1.162+.01 

.0646+. 007 

9.0 

1.00009 

.150X10"5 

4.0 

1.093+.01 

.0353+. 005 
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Table  15 

Pair  Distribution  Function:  Three-Particle  Pair  Exchange 
Contribution  for  T  =  1 . 7°K 


r 

U(r) 

n_    (r)/X3 
2 

r 

U(r) 

np    (r)/X3 
2 

0.9 

.31+. 04 

.0231+. 006 

2.40 

1. 554+. 04 

.0423+. 006 

1.0 

.63+.1 

.012+. 02 

2.556 

1.521+.03 

.0472+. 005 

1.10 

1.34+.1 

.051+. 02 

2.778 

1.391+.03 

.0405+. 005 

1.15 

1.45+.1 

.032+. 02 

3.00 

1. 283+. 01 
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